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;
}