Finding linear regression bearing from set of coordinates isn't giving desired result

I have a series of coordinates that I'm wanting to find the linear regression for, so I can then find the bearing of the line. I'm using the Linear Regression algorithm from Swift Algorithm Club on this set of coordinates:

51.48163827836369, -0.019464668521006683
51.481654860705667, -0.019452641350085287
51.481674140657908, -0.01943882290242982
51.481690344713748, -0.019423713183982727
51.481705506128442, -0.019419045258473489
51.481722553489625, -0.01940979681751287
51.48173752576799, -0.019412687104136239
51.48174912673462, -0.019409150632213823
51.4817646359283, -0.019389300889997685
51.481779676567427, -0.019388957697628939
51.481792568262044, -0.019402453532393338
51.481804168699682, -0.019415663863242116
51.481822746271966, -0.019423725406568337
51.481838162880258, -0.019428618620622728
51.481855587496689, -0.01942372804705883
51.481867836051975, -0.019430554178484272
51.481883136496599, -0.019432972610502475
51.481899505688553, -0.019425501321734373
51.481914919015246, -0.019424832166464512
51.481932613348015, -0.019457392982985766
51.481949346615657, -0.019472056279255412
51.481968002664601, -0.019458232757642691
51.481986902059589, -0.019446792660346546
51.482003086257393, -0.019433403642779012

Adding this to a map using darrinward.com, I get this path, moving roughly north:

So from there, to get the bearing I take the regression of the first coordinate's longitude (to get the latitude), :

    let regression = linearRegression(longitudes, latitudes)

    let firstRegressionCoordinate = CLLocationCoordinate2D(latitude: regression(firstCoordinate.longitude), longitude:firstCoordinate.longitude)
    let lastRegressionCoordinate = CLLocationCoordinate2D(latitude: regression(lastCoordinate.longitude), longitude: lastCoordinate.longitude)

    return firstRegressionCoordinate.bearing(to: lastRegressionCoordinate)

My bearing function works properly, and gives me a bearing of 156.20969º. What I'm expecting is a number closer to 0º, or closer to 360º, given the actual path taken.

I feel the issue isn't with the linear regression algorithm, or with the bearing function, but with my coordinates I'm asking for in order to determine the bearing. Am I asking it for the wrong thing? Is there something else I should be doing differently?

The complete code, that can be run in a playground:

import CoreLocation

public extension FloatingPoint {
    public var degreesToRadians: Self { return self * .pi / 180 }
    public var radiansToDegrees: Self { return self * 180 / .pi }
}

extension CLLocationCoordinate2D {
    ///Returns the initial bearing of travel to another coordinate
    func bearing(to: CLLocationCoordinate2D) -> CLLocationDegrees {
        let fromLatRadians = latitude.degreesToRadians
        let fromLongRadians = longitude.degreesToRadians

        let toLatRadians = to.latitude.degreesToRadians
        let toLongRadians = to.longitude.degreesToRadians

        let y = sin(toLongRadians - fromLongRadians) * cos(toLatRadians)
        let x = cos(fromLatRadians) * sin(toLatRadians) - sin(fromLatRadians) * cos(toLatRadians) * cos(toLongRadians - fromLongRadians)

        var bearing = atan2(y, x).radiansToDegrees

        bearing = (bearing + 360.0).truncatingRemainder(dividingBy: 360.0)

        return bearing
    }
}

extension Array where Element == CLLocationCoordinate2D {
    func linearRegressionBearing() -> Double {
        var longitudes = [CLLocationDegrees]()
        var latitudes = [CLLocationDegrees]()

        for coordinate in self {
            longitudes.append(coordinate.longitude)
            latitudes.append(coordinate.latitude)
        }

        let regression = linearRegression(longitudes, latitudes)

        let firstCoordinate = CLLocationCoordinate2D(latitude: regression(self.first!.longitude), longitude: self.first!.longitude)
        let lastCoordinate = CLLocationCoordinate2D(latitude: regression(self.last!.longitude), longitude: self.last!.longitude)

        return firstCoordinate.bearing(to: lastCoordinate)
    }
}

// A closed form solution
func average(_ input: [Double]) -> Double {
    return input.reduce(0, +) / Double(input.count)
}

func multiply(_ a: [Double], _ b: [Double]) -> [Double] {
    return zip(a, b).map(*)
}

func linearRegression(_ xs: [Double], _ ys: [Double]) -> (Double) -> Double {
    let sum1 = average(multiply(xs, ys)) - average(xs) * average(ys)
    let sum2 = average(multiply(xs, xs)) - pow(average(xs), 2)
    let slope = sum1 / sum2
    let intercept = average(ys) - slope * average(xs)
    return { x in intercept + slope * x }
}

let coordinates = [
    CLLocationCoordinate2D(latitude: 51.48163827836369, longitude: -0.019464668521006683),
    CLLocationCoordinate2D(latitude: 51.481654860705667, longitude: -0.019452641350085287),
    CLLocationCoordinate2D(latitude: 51.481674140657908, longitude: -0.01943882290242982),
    CLLocationCoordinate2D(latitude: 51.481690344713748, longitude: -0.019423713183982727),
    CLLocationCoordinate2D(latitude: 51.481705506128442, longitude: -0.019419045258473489),
    CLLocationCoordinate2D(latitude: 51.481722553489625, longitude: -0.01940979681751287),
    CLLocationCoordinate2D(latitude: 51.48173752576799, longitude: -0.019412687104136239),
    CLLocationCoordinate2D(latitude: 51.48174912673462, longitude: -0.019409150632213823),
    CLLocationCoordinate2D(latitude: 51.4817646359283, longitude: -0.019389300889997685),
    CLLocationCoordinate2D(latitude: 51.481779676567427, longitude: -0.019388957697628939),
    CLLocationCoordinate2D(latitude: 51.481792568262044, longitude: -0.019402453532393338),
    CLLocationCoordinate2D(latitude: 51.481804168699682, longitude: -0.019415663863242116),
    CLLocationCoordinate2D(latitude: 51.481822746271966, longitude: -0.019423725406568337),
    CLLocationCoordinate2D(latitude: 51.481838162880258, longitude: -0.019428618620622728),
    CLLocationCoordinate2D(latitude: 51.481855587496689, longitude: -0.01942372804705883),
    CLLocationCoordinate2D(latitude: 51.481867836051975, longitude: -0.019430554178484272),
    CLLocationCoordinate2D(latitude: 51.481883136496599, longitude: -0.019432972610502475),
    CLLocationCoordinate2D(latitude: 51.481899505688553, longitude: -0.019425501321734373),
    CLLocationCoordinate2D(latitude: 51.481914919015246, longitude: -0.019424832166464512),
    CLLocationCoordinate2D(latitude: 51.481932613348015, longitude: -0.019457392982985766),
    CLLocationCoordinate2D(latitude: 51.481949346615657, longitude: -0.019472056279255412),
    CLLocationCoordinate2D(latitude: 51.481968002664601, longitude: -0.019458232757642691),
    CLLocationCoordinate2D(latitude: 51.481986902059589, longitude: -0.019446792660346546),
    CLLocationCoordinate2D(latitude: 51.482003086257393, longitude: -0.019433403642779012)
]

coordinates.linearRegressionBearing()

I confirm that your bearing function works correctly. However, the regression value you get for latitude is 51.45437262420372 - which is where the regression line intersects with longitude = 0.0. I think you need to use both the intersection point and slope of the regression line to calculate fitted lat,lon values (within your map region), and then find the bearing using those values.

LATER ADDITION: An important issue here is that longitude coordinates are not linear in relation to distance - the distance covered by a degree of longitude varies with latitude. I think the simplest solution is just to convert your array of lat,lon values to lat,normalisedLon values where normalisedLon = lon * cosine(lat).

When you do this the angle of your path from your first coordinate to your last coordinate = 3.055 degrees (just east of north which appears to be correct).

Once you've done this then the slope of your regression line will in fact represent the direction you are seeking i.e. heading = atan(slope).

However, your regression code just doesn't seem to work for me. It should give a result very close to 3 degrees also, but it is way off.