Failing to calculate the gradient and y-shift for best fit line using sum-squares method. Can anyone help me to pinpoint the issue?

Question

I'm trying to work out how to make the working implementation of fitting the straight line among the set data points. I base my function:

pub fn sum_squares<Q: CanBeLabel>(l1: &Line<Q>, accuracy: u32) -> Point {

    let mut dataset_size: f32 = 0.;
    let (mut y_mean, mut x_mean): (f32, f32) = (0., 0.);
    let (mut y_max, mut x_max): (f32, f32) = (0., 0.);

    for point in l1.points() {
        dataset_size += 1.;
        y_mean += point.y;
        x_mean += point.x;

        if y_max < point.y {
            y_max = point.y;
        }

        if x_max < point.x {
            x_max = point.x;
        }
    }

    y_mean /= dataset_size;
    x_mean /= dataset_size;

    let fpa = |a: f32, b: f32| -> f32 {
        2. * dataset_size * a * x_mean * x_mean + 2. * dataset_size * x_mean * b - 2. * dataset_size * y_mean * x_mean
    };

    let fpb = |a: f32, b: f32| -> f32 {
        - (a * dataset_size * y_mean) + 2. * a * k * x_mean + 2. * dataset_size * b
    };

    let fpa2 = |b: f32| -> f32 { 2. * dataset_size * x_mean  };

    let fpb2 = |a: f32| -> f32 { 2. * dataset_size  };

    let mut a: f32 = y_max;
    let mut b: f32 = y_max;

    for _ in 0..accuracy {
        a -=  ( fpa(a, b) / fpa2(b));
        b -= (fpb(a, b) / fpb2(a));

    }

    return Point::new(a, b);
}

I base my function on the following expansion:

$\sum^{k}_{n = 1}\left(y_n-\left(ax_n+b\right)\right)^2=$

$\sum^{k}_{n = 1}{{(y}_n^2-2y_n\left(ax_n+b\right)+\left(ax_n+b\right)^2)}=$

$\sum^{k}_{n = 1}\left(y_n^2\ -2y_nax_n-2y_nb+a^2{x_n}^2+2a^2x_nb+b^2\right)=$

$k\bar{y^2}\ -2k\bar{y}ax-2k\bar{y}b+{a^2k\bar{x}}^2+2a^2k\bar{x}b+kb^2$

and then take derivatives with respect to a, and b:

a. $0\ -\ 2k\bar{y}x\ -\ 0+{\ 2ak\bar{x}}^2\ +\ 4ak\bar{x}b\ \ +\ 0\ \equiv{\ 2ak\bar{x}}^2\ +\ \ 4ak\bar{x}b\ -\ 2k\bar{y}x$

b. $0\ -\ 0\ -\ 2k\bar{y}b\ +\ 0\ +\ 2a^2k\bar{x}b\ +\ kb^2\ \equiv2a^2k\bar{x}\ +\ 2kb\ -\ 2k\bar{y}b$

and then approximating the zeroes of the derivative via Newton-Raphson (so I also use second derivatives, where $2k\bar{x}^2 + 4k\bar{x}b$ is with respect to a and $2k - 2k\bar{y}$ with respect to b).

However when I plot my graph for a dataset consisting of 4 points where $x \in [1,4] \land x \in \mathbb{Z}$ and $y = x + 666$,

fn make_dataset() -> Line<String> {
    let mut ret: Line<String> = Line::new("Cats".to_owned());

    for i in 0..5 {
        ret.add_point(Point::new(i as f32, i as f32 + 7_f32));
    }

    ret
}

and for checking:

def make_line(a: float, b: float, xs: range):
    x, y = [], []
    for i in xs:
        x.append(i), y.append(a + 666)

    return x, y

for every pair of consequetive x values graph looks like so:

for odd x's:

for even x's:

Answer 1

Just like @gnasher729, I did not check your code nor your equations, which seem overly complicated. An explicit solution is

$$a=\frac{\sum_k (y_k-\overline y)(x_k-\overline x)}{\sum_k (x-\overline x)^2}$$

and

$$b=\overline y-a\overline x.$$

---

Explanation:

Note that by linearity (or by differentiating the global error on $b$),

$$\overline y=a\overline x+b$$ allows to eliminate $b$ and reduce the problem to the minimization of

$$\sum_k(\hat y_k-a\hat x_k)^2$$ where the hat denotes centered variables. After differentiation on $a$,

$$\sum_k\hat y_k\hat x_k-a\sum_k\hat x_k^2=0.$$

Answer 2

I really don't want to go and debug your code. But if you want a linear approximation by arbitrary functions $f_i$ minimising the squares of errors, you create a system of linear equations, with the coefficients on the left side being the sums of $f_i \cdot f_j$, and on the right side you write the sums of $f_i \cdot y$. You solve the system and that's it.

Accordingly, for approximation by two functions $f_1(x) = x$ and $f_2(x) = 1$ you create a system of two linear equations in two variables, with the four coefficients on the left side being the sums of x^2, x, x and 1, and the numbers on the right side being the sums of $x \cdot y$ and $1 \cdot y$.

This is all just linear algebra and linear equations, so using Newton-Raphson seems quite bizarre.