The best fit bearing should be sufficiently approximated as the average azimuth over the sequence:
SELECT AVG(azm) AS bearing
SELECT degrees(ST_Azimuth(geom, LEAD(geom) OVER(ORDER BY <sequence_id>))) AS azm
Alternatively, you can easily reproduce the approximation as described in the linked answer with PostgreSQL's
SELECT ATAN2D(1, (REGR_SLOPE(ST_Y(geom), ST_X(geom)) + REGR_SLOPE(ST_X(geom), ST_Y(geom)) / 2) AS bearing
(Switched ordinates in
ATAN2D to get angle relative to true north)
The regression analysis is stronger influenced by outliers, the averaging tends to have steeper angles.
From a strict statistical POV none of these are the best-fit line as per Deming; the regression line average is closer, but completely unaware of the spheroidal character of the 'measures' (thus more precise on a 2D cartesian plane).