<![CDATA[
#include <stdio.h>
#include <numeric.h>
#define NPT 100 /* Number of data points */
#define NTERM 5 /* Number sof term */
#define SPREAD 0.1 /* Number sof term */
int main() {
int i,j,status,ia[NTERM];
array double u[NPT],x[NPT],y[NPT],a[NTERM],dev[NTERM];
array double sig[NPT],covar[NTERM][NTERM];
double chisq;
/* Create a data set of NTERM order polynomial with gaussian deviation*/
linspace(x,0.1,0.1*NPT);
y = 8*x.*x.*x.*x + 5*x.*x.*x + 3*x.*x + 6*x + (array double[NPT])(7);
/* put uniform random deviation on data set */
urand(u);
y += 0.1*u;
sig=(array double[NPT])(SPREAD);
ia[0] = 1; ia[1] = 0; ia[2] = 1; ia[3] = 0, ia[4] = 1;
a[1] = 5; a[3] = 6;
status=polyfit(a,x,y,sig,ia,covar,&chisq);
dev = sqrt(diagonal(covar));
if(status) printf("Abnormal fit");
printf("\n%11s%23s\n","Coefficients");
printf("\n%11s %21s\n","parameter","uncertainty");
for (i=0;i<NTERM;i++)
printf(" a[%1d] = %8.6f %12.6f\n",
i,a[i],dev[i]);
printf("chi-square = %12f\n",chisq);
printf("full covariance matrix\n");
printf("%12f",covar);
}
]]>
<![CDATA[
Coefficients
parameter uncertainty
a[0] = 8.000007 0.000013
a[1] = 5.000000 0.000000
a[2] = 2.999413 0.001165
a[3] = 6.000000 0.000000
a[4] = 7.057378 0.018870
chi-square = 7.414093
full covariance matrix
0.000000 0.000000 -0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000
-0.000000 0.000000 0.000001 0.000000 -0.000016
0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 -0.000016 0.000000 0.000356
]]>