- /*
- double sxy,sxx,syy,sx,sy,w,s_2,dx2,dy2,mins;
-
- sxy=sxx=syy=sx=sy=w=0.0;
- mins = dy[0];
- for(i=1; (i<n); i++)
- mins = min(mins,dy[i]);
- if (mins <= 0)
- gmx_fatal(FARGS,"Zero or negative weigths in linear regression analysis");
-
- for (i=0; i<n; i++) {
- s_2 = dsqr(1.0/dy[i]);
- sxx += s_2*dsqr(x[i]);
- sxy += s_2*y[i]*x[i];
- syy += s_2*dsqr(y[i]);
- sx += s_2*x[i];
- sy += s_2*y[i];
- w += s_2;
- }
- sxx = sxx/w;
- sxy = sxy/w;
- syy = syy/w;
- sx = sx/w;
- sy = sy/w;
- dx2 = (sxx-sx*sx);
- dy2 = (syy-sy*sy);
- * a=(sxy-sy*sx)/dx2;
- * b=(sy-(*a)*sx);
-
- * chi2=0;
- for(i=0; i<n; i++)
- * chi2+=dsqr((y[i]-((*a)*x[i]+(*b)))/dy[i]);
- * chi2 = *chi2/w;
-
- * da = sqrt(*chi2/((n-2)*dx2));
- * db = *da*sqrt(sxx);
- * r = *a*sqrt(dx2/dy2);
-
- if (debug)
- fprintf(debug,"sx = %g, sy = %g, sxy = %g, sxx = %g, w = %g\n"
- "chi2 = %g, dx2 = %g\n",
- sx,sy,sxy,sxx,w,*chi2,dx2);
-
- if (n > 2)
- * chi2 = sqrt(*chi2/(n-2));
- else
- * chi2 = 0;
- */