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