A program to fit straight lines to data

#include <math.h>
#include <stdio.h>
#include <assert.h>

/*
Michael Ashley / UNSW / 06-Jun-2004
*/

void linfit(double x[], double y[], double sy[], int n,
double *a, double *b, double *sa, double *sb) {

/*
Fits a straight line y = a + bx to the n data pairs (x,y).
sy is the standard deviation of y. The standard deviations
of a and b are also returned.
*/

double S=0.0, Sx=0.0, Sy=0.0, Stt=0.0, Sty=0.0;
double t;
int i;

assert(n > 1);

for (i = 0; i < n; i++) {
double ss;
ss = sy[i]*sy[i];
S += 1/ss;
Sx += x[i]/ss;
Sy += y[i]/ss;
}

for (i = 0; i < n; i++) {
t = (x[i]-Sx/S)/sy[i];
Stt += t*t;
Sty += t*y[i]/sy[i];
}

*b = Sty / Stt;
*a = (Sy - Sx* *b)/S;
*sa = sqrt((1+(Sx*Sx)/(S*Stt))/S);
*sb = sqrt(1/Stt);
}

int main(void) {
double x[3]={1.0, 2.0, 3.0};
double y[3]={1.1, 2.3, 3.5};
double sy[3]={0.01, 0.01, 0.01};

double a, b, sa, sb;

linfit(x, y, sy, (sizeof x)/(sizeof x[0]), &a, &b, &sa, &sb);

printf("y = (%lf+/-%lf) + (%lf+/-%lf) * x\n", a, sa, b, sb);
return 0;
}