src/tools/expfit.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42
  43 #include <sysstuff.h>
  44 #include <string.h>
  45 #include <math.h>
  46 #include "typedefs.h"
  47 #include "smalloc.h"
  48 #include "xvgr.h"
  49 #include "futil.h"
  50 #include "gstat.h"
  51 #include "vec.h"
  52 #include "statutil.h"
  53 #include "index.h"
  54
  55 const int  nfp_ffn[effnNR] = { 0, 1, 2, 3, 2, 5, 7, 9, 4, 3};
  56
  57 const char *s_ffn[effnNR+2] = { NULL, "none", "exp", "aexp", "exp_exp", "vac",
  58                           "exp5", "exp7", "exp9","erffit", NULL, NULL };
  59 /* We don't allow errest as a choice on the command line */
  60
  61 const char *longs_ffn[effnNR] = {
  62   "no fit",
  63   "y = exp(-x/a1)",
  64   "y = a2 exp(-x/a1)",
  65   "y = a2 exp(-x/a1) + (1-a2) exp(-x/a3)",
  66   "y = exp(-v) (cosh(wv) + 1/w sinh(wv)), v = x/(2 a1), w = sqrt(1 - a2)",
  67   "y = a1 exp(-x/a2) +  a3 exp(-x/a4) + a5",
  68   "y = a1 exp(-x/a2) +  a3 exp(-x/a4) + a5 exp(-x/a6) + a7",
  69   "y = a1 exp(-x/a2) +  a3 exp(-x/a4) + a5 exp(-x/a6) + a7 exp(-x/a8) + a9",
  70   "y = 1/2*(a1+a2) - 1/2*(a1-a2)*erf( (x-a3) /a4^2)",
  71   "y = a2*ee(a1,x) + (1-a2)*ee(a2,x)"
  72 };
  73
  74 extern gmx_bool mrqmin_new(real x[],real y[],real sig[],int ndata,real a[],
  75                        int ia[],int ma,real **covar,real **alpha,real *chisq,
  76                        void (*funcs)(real, real [], real *, real []),
  77                        real *alamda);
  78
  79 static real myexp(real x,real A,real tau)
  80 {
  81   if ((A == 0) || (tau == 0))
  82     return 0;
  83   return A*exp(-x/tau);
  84 }
  85
  86 void erffit (real x, real a[], real *y, real dyda[])
  87 {
  88 /* Fuction
  89  *      y=(a1+a2)/2-(a1-a2)/2*erf((x-a3)/a4^2)
  90  */
  91
  92 real erfarg;
  93 real erfval;
  94 real erfarg2;
  95 real derf;
  96
  97  erfarg=(x-a[3])/(a[4]*a[4]);
  98         erfarg2=erfarg*erfarg;
  99         erfval=gmx_erf(erfarg)/2;
 100         derf=M_2_SQRTPI*(a[1]-a[2])/2*exp(-erfarg2)/(a[4]*a[4]);
 101         *y=(a[1]+a[2])/2-(a[1]-a[2])*erfval;
 102         dyda[1]=1/2-erfval;
 103         dyda[2]=1/2+erfval;
 104         dyda[3]=derf;
 105         dyda[4]=2*derf*erfarg;
 106 }
 107
 108
 109
 110 static void exp_one_parm(real x,real a[],real *y,real dyda[])
 111 {
 112   /* Fit to function
 113    *
 114    * y = exp(-x/a1)
 115    *
 116    */
 117
 118   real e1;
 119
 120   e1      = exp(-x/a[1]);
 121   *y      = e1;
 122   dyda[1] = x*e1/sqr(a[1]);
 123 }
 124
 125 static void exp_two_parm(real x,real a[],real *y,real dyda[])
 126 {
 127   /* Fit to function
 128    *
 129    * y = a2 exp(-x/a1)
 130    *
 131    */
 132
 133   real e1;
 134
 135   e1      = exp(-x/a[1]);
 136   *y      = a[2]*e1;
 137   dyda[1] = x*a[2]*e1/sqr(a[1]);
 138   dyda[2] = e1;
 139 }
 140
 141 static void exp_3_parm(real x,real a[],real *y,real dyda[])
 142 {
 143   /* Fit to function
 144    *
 145    * y = a2 exp(-x/a1) + (1-a2) exp(-x/a3)
 146    *
 147    */
 148
 149   real e1,e2;
 150
 151   e1      = exp(-x/a[1]);
 152   e2      = exp(-x/a[3]);
 153   *y      = a[2]*e1 + (1-a[2])*e2;
 154   dyda[1] = x*a[2]*e1/sqr(a[1]);
 155   dyda[2] = e1-e2;
 156   dyda[3] = x*(1-a[2])*e2/sqr(a[3]);
 157 }
 158
 159 static void exp_5_parm(real x,real a[],real *y,real dyda[])
 160 {
 161   /* Fit to function
 162    *
 163    * y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5
 164    *
 165    */
 166
 167   real e1,e2;
 168
 169   e1      = exp(-x/a[2]);
 170   e2      = exp(-x/a[4]);
 171   *y      = a[1]*e1 + a[3]*e2 + a[5];
 172
 173   if (debug)
 174     fprintf(debug,"exp_5_parm called: x = %10.3f  y = %10.3f\n"
 175             "a = ( %8.3f  %8.3f  %8.3f  %8.3f  %8.3f)\n",
 176             x,*y,a[1],a[2],a[3],a[4],a[5]);
 177   dyda[1] = e1;
 178   dyda[2] = x*e1/sqr(a[2]);
 179   dyda[3] = e2;
 180   dyda[4] = x*e2/sqr(a[4]);
 181   dyda[5] = 0;
 182 }
 183
 184 static void exp_7_parm(real x,real a[],real *y,real dyda[])
 185 {
 186   /* Fit to function
 187    *
 188    * y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7
 189    *
 190    */
 191
 192   real e1,e2,e3;
 193
 194   e1      = exp(-x/a[2]);
 195   e2      = exp(-x/a[4]);
 196   e3      = exp(-x/a[6]);
 197   *y      = a[1]*e1 + a[3]*e2 + a[5]*e3 + a[7];
 198
 199   dyda[1] = e1;
 200   dyda[2] = x*e1/sqr(a[2]);
 201   dyda[3] = e2;
 202   dyda[4] = x*e2/sqr(a[4]);
 203   dyda[5] = e3;
 204   dyda[6] = x*e3/sqr(a[6]);
 205   dyda[7] = 0;
 206 }
 207
 208 static void exp_9_parm(real x,real a[],real *y,real dyda[])
 209 {
 210   /* Fit to function
 211    *
 212    * y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7
 213    *
 214    */
 215
 216   real e1,e2,e3,e4;
 217
 218   e1      = exp(-x/a[2]);
 219   e2      = exp(-x/a[4]);
 220   e3      = exp(-x/a[6]);
 221   e4      = exp(-x/a[8]);
 222   *y      = a[1]*e1 + a[3]*e2 + a[5]*e3 + a[7]*e4 + a[9];
 223
 224   dyda[1] = e1;
 225   dyda[2] = x*e1/sqr(a[2]);
 226   dyda[3] = e2;
 227   dyda[4] = x*e2/sqr(a[4]);
 228   dyda[5] = e3;
 229   dyda[6] = x*e3/sqr(a[6]);
 230   dyda[7] = e4;
 231   dyda[8] = x*e4/sqr(a[8]);
 232   dyda[9] = 0;
 233 }
 234
 235 static void vac_2_parm(real x,real a[],real *y,real dyda[])
 236 {
 237   /* Fit to function
 238    *
 239    * y = 1/2 (1 - 1/w) exp(-(1+w)v) + 1/2 (1 + 1/w) exp(-(1-w)v)
 240    *
 241    *   = exp(-v) (cosh(wv) + 1/w sinh(wv))
 242    *
 243    *    v = x/(2 a1)
 244    *    w = sqrt(1 - a2)
 245    *
 246    *    For tranverse current autocorrelation functions:
 247    *       a1 = tau
 248    *       a2 = 4 tau (eta/rho) k^2
 249    *
 250    */
 251
 252   double v,det,omega,omega2,em,ec,es;
 253
 254   v   = x/(2*a[1]);
 255   det = 1 - a[2];
 256   em  = exp(-v);
 257   if (det != 0) {
 258     omega2  = fabs(det);
 259     omega   = sqrt(omega2);
 260     if (det > 0) {
 261       ec = em*0.5*(exp(omega*v)+exp(-omega*v));
 262       es = em*0.5*(exp(omega*v)-exp(-omega*v))/omega;
 263     } else {
 264       ec = em*cos(omega*v);
 265       es = em*sin(omega*v)/omega;
 266     }
 267     *y      = ec + es;
 268     dyda[2] = (v/det*ec+(v-1/det)*es)/(-2.0);
 269     dyda[1] = (1-det)*v/a[1]*es;
 270   } else {
 271     *y      = (1+v)*em;
 272     dyda[2] = -v*v*em*(0.5+v/6);
 273     dyda[1] = v*v/a[1]*em;
 274   }
 275 }
 276
 277 static void errest_3_parm(real x,real a[],real *y,real dyda[])
 278 {
 279   real e1,e2,v1,v2;
 280
 281   if (a[1])
 282     e1 = exp(-x/a[1]) - 1;
 283   else
 284     e1 = 0;
 285   if (a[3])
 286     e2 = exp(-x/a[3]) - 1;
 287   else
 288     e2 = 0;
 289
 290   if (x > 0) {
 291     v1 = 2*a[1]*(e1*a[1]/x + 1);
 292     v2 = 2*a[3]*(e2*a[3]/x + 1);
 293     *y      = a[2]*v1 + (1-a[2])*v2;
 294     dyda[1] = 2*     a[2] *(v1/a[1] + e1);
 295     dyda[3] = 2*(1 - a[2])*(v2/a[3] + e2);
 296     dyda[2] = (v1 - v2);
 297   } else {
 298     *y      = 0;
 299     dyda[1] = 0;
 300     dyda[3] = 0;
 301     dyda[2] = 0;
 302   }
 303 }
 304
 305 typedef void (*myfitfn)(real x,real a[],real *y,real dyda[]);
 306 myfitfn mfitfn[effnNR] = {
 307   exp_one_parm, exp_one_parm, exp_two_parm, exp_3_parm, vac_2_parm,
 308   exp_5_parm,   exp_7_parm,   exp_9_parm, erffit,  errest_3_parm
 309 };
 310
 311 real fit_function(int eFitFn,real *parm,real x)
 312 {
 313   static real y,dum[8];
 314
 315   mfitfn[eFitFn](x,parm-1,&y,dum);
 316
 317   return y;
 318 }
 319
 320 /* lmfit_exp supports up to 3 parameter fitting of exponential functions */
 321 static gmx_bool lmfit_exp(int nfit,real x[],real y[],real dy[],real ftol,
 322                       real parm[],real dparm[],gmx_bool bVerbose,
 323                       int eFitFn,int fix)
 324 {
 325   real chisq,ochisq,alamda;
 326   real *a,**covar,**alpha,*dum;
 327   gmx_bool bCont;
 328   int  i,j,ma,mfit,*lista,*ia;
 329
 330   if ((eFitFn < 0) || (eFitFn >= effnNR))
 331     gmx_fatal(FARGS,"fitfn = %d, should be in 0..%d (%s,%d)",
 332                 effnNR-1,eFitFn,__FILE__,__LINE__);
 333
 334   ma=mfit=nfp_ffn[eFitFn];         /* number of parameters to fit */
 335   snew(a,ma+1);
 336   snew(covar,ma+1);
 337   snew(alpha,ma+1);
 338   snew(lista,ma+1);
 339   snew(ia,ma+1);
 340   snew(dum,ma+1);
 341   for(i=1; (i<ma+1); i++) {
 342     lista[i] = i;
 343     ia[i] = 1; /* fixed bug B.S.S 19/11  */
 344     snew(covar[i],ma+1);
 345     snew(alpha[i],ma+1);
 346   }
 347   if (fix) {
 348     if (bVerbose)
 349       fprintf(stderr,"Will keep parameters fixed during fit procedure: %d\n",
 350               fix);
 351     for(i=0; i<ma; i++)
 352       if (fix & 1<<i)
 353         ia[i+1] = 0;
 354   }
 355   if (debug)
 356     fprintf(debug,"%d parameter fit\n",mfit);
 357
 358   /* Initial params */
 359   alamda = -1;    /* Starting value   */
 360   chisq  = 1e12;
 361   for(i=0; (i<mfit); i++)
 362     a[i+1] = parm[i];
 363
 364   j = 0;
 365   if (bVerbose)
 366     fprintf(stderr,"%4s  %10s  %10s  %10s  %10s  %10s %10s\n",
 367             "Step","chi^2","Lambda","A1","A2","A3","A4");
 368   do {
 369     ochisq = chisq;
 370     /* mrqmin(x-1,y-1,dy-1,nfit,a,ma,lista,mfit,covar,alpha,
 371      *   &chisq,expfn[mfit-1],&alamda)
 372      */
 373     if (!mrqmin_new(x-1,y-1,dy-1,nfit,a,ia,ma,covar,alpha,&chisq,
 374                     mfitfn[eFitFn],&alamda))
 375       return FALSE;
 376
 377     if (bVerbose) {
 378       fprintf(stderr,"%4d  %10.5e  %10.5e  %10.5e",
 379               j,chisq,alamda,a[1]);
 380       if (mfit > 1)
 381         fprintf(stderr,"  %10.5e",a[2]);
 382       if (mfit > 2)
 383         fprintf(stderr,"  %10.5e",a[3]);
 384       if (mfit > 3)
 385          fprintf(stderr," %10.5e",a[4]);
 386       fprintf(stderr,"\n");
 387     }
 388     j++;
 389     bCont = ((fabs(ochisq - chisq) > fabs(ftol*chisq)) ||
 390              ((ochisq == chisq)));
 391   } while (bCont && (alamda != 0.0) && (j < 50));
 392   if (bVerbose)
 393     fprintf(stderr,"\n");
 394
 395   /* Now get the covariance matrix out */
 396   alamda = 0;
 397
 398   /*  mrqmin(x-1,y-1,dy-1,nfit,a,ma,lista,mfit,covar,alpha,
 399    * &chisq,expfn[mfit-1],&alamda)
 400    */
 401   if ( !mrqmin_new(x-1,y-1,dy-1,nfit,a,ia,ma,covar,alpha,&chisq,
 402                    mfitfn[eFitFn],&alamda))
 403     return FALSE;
 404
 405   for(j=0; (j<mfit); j++) {
 406     parm[j]  = a[j+1];
 407     dparm[j] = covar[j+1][j+1];
 408   }
 409
 410   for(i=0; (i<ma+1); i++) {
 411     sfree(covar[i]);
 412     sfree(alpha[i]);
 413   }
 414   sfree(a);
 415   sfree(covar);
 416   sfree(alpha);
 417   sfree(lista);
 418   sfree(dum);
 419
 420   return TRUE;
 421 }
 422
 423 real do_lmfit(int ndata,real c1[],real sig[],real dt,real x0[],
 424               real begintimefit,real endtimefit,const output_env_t oenv,
 425               gmx_bool bVerbose, int eFitFn,real fitparms[],int fix)
 426 {
 427   FILE *fp;
 428   char buf[32];
 429
 430   int  i,j,nparm,nfitpnts;
 431   real integral,ttt;
 432   real *parm,*dparm;
 433   real *x,*y,*dy;
 434   real ftol = 1e-4;
 435
 436   nparm = nfp_ffn[eFitFn];
 437   if (debug) {
 438     fprintf(debug,"There are %d points to fit %d vars!\n",ndata,nparm);
 439     fprintf(debug,"Fit to function %d from %g through %g, dt=%g\n",
 440             eFitFn,begintimefit,endtimefit,dt);
 441   }
 442
 443   snew(x,ndata);
 444   snew(y,ndata);
 445   snew(dy,ndata);
 446
 447   j=0;
 448   for(i=0; i<ndata; i++) {
 449     ttt = x0 ? x0[i] : dt*i;
 450     if (ttt>=begintimefit && ttt<=endtimefit) {
 451       x[j] = ttt;
 452       y[j] = c1[i];
 453
 454       /* mrqmin does not like sig to be zero */
 455       if (sig[i]<1.0e-7)
 456         dy[j]=1.0e-7;
 457       else
 458         dy[j]=sig[i];
 459       if (debug)
 460         fprintf(debug,"j= %d, i= %d, x= %g, y= %g, dy= %g\n",
 461                 j,i,x[j],y[j],dy[j]);
 462       j++;
 463     }
 464   }
 465   nfitpnts = j;
 466   integral = 0;
 467   if (nfitpnts < nparm)
 468     fprintf(stderr,"Not enough data points for fitting!\n");
 469   else {
 470     snew(parm,nparm);
 471     snew(dparm,nparm);
 472     if (fitparms)
 473       for(i=0; (i < nparm); i++)
 474         parm[i]=fitparms[i];
 475
 476     if (!lmfit_exp(nfitpnts,x,y,dy,ftol,parm,dparm,bVerbose,eFitFn,fix))
 477       fprintf(stderr,"Fit failed!\n");
 478     else if (nparm <= 3) {
 479       /* Compute the integral from begintimefit */
 480       if (nparm == 3)
 481         integral=(parm[0]*myexp(begintimefit,parm[1],  parm[0]) +
 482                   parm[2]*myexp(begintimefit,1-parm[1],parm[2]));
 483       else if (nparm == 2)
 484         integral=parm[0]*myexp(begintimefit,parm[1],  parm[0]);
 485       else if (nparm == 1)
 486         integral=parm[0]*myexp(begintimefit,1,  parm[0]);
 487       else
 488         gmx_fatal(FARGS,"nparm = %d in file %s, line %d",
 489                     nparm,__FILE__,__LINE__);
 490
 491       /* Generate THE output */
 492       if (bVerbose) {
 493         fprintf(stderr,"FIT: # points used in fit is: %d\n",nfitpnts);
 494         fprintf(stderr,"FIT: %21s%21s%21s\n",
 495                 "parm0     ","parm1 (ps)   ","parm2 (ps)    ");
 496         fprintf(stderr,"FIT: ------------------------------------------------------------\n");
 497         fprintf(stderr,"FIT: %8.3g +/- %8.3g%9.4g +/- %8.3g%8.3g +/- %8.3g\n",
 498                 parm[0],dparm[0],parm[1],dparm[1],parm[2],dparm[2]);
 499         fprintf(stderr,"FIT: Integral (calc with fitted function) from %g ps to inf. is: %g\n",
 500                 begintimefit,integral);
 501
 502         sprintf(buf,"test%d.xvg",nfitpnts);
 503         fp = xvgropen(buf,"C(t) + Fit to C(t)","Time (ps)","C(t)",oenv);
 504         fprintf(fp,"# parm0 = %g, parm1 = %g, parm2 = %g\n",
 505                 parm[0],parm[1],parm[2]);
 506         for(j=0; (j<nfitpnts); j++) {
 507           ttt = x0 ? x0[j] : dt*j;
 508           fprintf(fp,"%10.5e  %10.5e  %10.5e\n",
 509                   ttt,c1[j],fit_function(eFitFn,parm,ttt));
 510         }
 511         xvgrclose(fp);
 512       }
 513     }
 514     if (fitparms)
 515     {
 516         for(i=0;(i<nparm);i++)
 517         {
 518             fitparms[i] = parm[i];
 519         }
 520     }
 521     sfree(parm);
 522     sfree(dparm);
 523   }
 524
 525   sfree(x);
 526   sfree(y);
 527   sfree(dy);
 528
 529   return integral;
 530 }
 531
 532 void do_expfit(int ndata,real c1[],real dt,real begintimefit,real endtimefit)
 533 {
 534   int i,n;
 535   real *x,*y,*Dy;
 536   real aa,bb,saa,sbb,A,tau,dA,dtau;
 537
 538   fprintf(stderr,"Will fit data from %g (ps) to %g (ps).\n",
 539           begintimefit,endtimefit);
 540
 541   snew(x,ndata);   /* allocate the maximum necessary space */
 542   snew(y,ndata);
 543   snew(Dy,ndata);
 544   n=0;
 545
 546   for(i=0; (i<ndata); i++) {
 547     if ( (dt*i >= begintimefit) && (dt*i <= endtimefit) ) {
 548       x[n]=dt*i;
 549       y[n]=c1[i];
 550       Dy[n]=0.5;
 551       fprintf(stderr,"n= %d, i= %d, x= %g, y= %g\n",n,i,x[n],y[n]);
 552       n++;
 553     }
 554   }
 555   fprintf(stderr,"# of data points used in the fit is : %d\n\n",n);
 556   expfit(n,x,y,Dy,&aa,&saa,&bb,&sbb);
 557
 558   A=exp(aa);
 559   dA=exp(aa)*saa;
 560   tau=-1.0/bb;
 561   dtau=sbb/sqr(bb);
 562   fprintf(stderr,"Fitted to y=exp(a+bx):\n");
 563   fprintf(stderr,"a = %10.5f\t b = %10.5f",aa,bb);
 564   fprintf(stderr,"\n");
 565   fprintf(stderr,"Fitted to y=Aexp(-x/tau):\n");
 566   fprintf(stderr,"A  = %10.5f\t tau  = %10.5f\n",A,tau);
 567   fprintf(stderr,"dA = %10.5f\t dtau = %10.5f\n",dA,dtau);
 568 }
 569
 570
 571 void expfit(int n, real *x, real *y, real *Dy, real *a, real *sa,
 572             real *b, real *sb)
 573 {
 574   real *w,*ly,A,SA,B,SB;
 575   int  i;
 576   real sum,xbar,ybar,Sxx,Sxy,wr2,chi2,gamma,Db;
 577
 578 #define ZERO 0.0
 579 #define ONE 1.0
 580 #define ONEP5 1.5
 581 #define TWO 2.0
 582
 583 #define sqr(x) ((x)*(x))
 584
 585   /*allocate memory */
 586   snew(w,n);
 587   snew(ly,n);
 588
 589   /* Calculate weights and values of ln(y). */
 590   for(i=0;(i<n); i++){
 591     w[i]=sqr(y[i]/Dy[i]);
 592     ly[i]=log(y[i]);
 593   }
 594
 595   /* The weighted averages of x and y: xbar and ybar. */
 596   sum=ZERO;
 597   xbar=ZERO;
 598   ybar=ZERO;
 599   for(i=0;(i<n);i++){
 600     sum+=w[i];
 601     xbar+=w[i]*x[i];
 602     ybar+=w[i]*ly[i];
 603   }
 604   xbar/=sum;
 605   ybar/=sum;
 606
 607   /* The centered product sums Sxx and Sxy, and hence A and B. */
 608   Sxx=ZERO;
 609   Sxy=ZERO;
 610   for(i=0;(i<n);i++){
 611     Sxx+=w[i]*sqr(x[i]-xbar);
 612     Sxy+=w[i]*(x[i]-xbar)*(ly[i]-ybar);
 613   }
 614   B=Sxy/Sxx;
 615   A=ybar-B*xbar;
 616
 617   /* Chi-squared (chi2) and gamma. */
 618   chi2=ZERO;
 619   gamma=ZERO;
 620   for(i=0;(i<n);i++){
 621     wr2=w[i]*sqr(ly[i]-A-B*x[i]);
 622     chi2+=wr2;
 623     gamma+=wr2*(x[i]-xbar);
 624   }
 625
 626   /* Refined values of A and B. Also SA and SB. */
 627   Db=-ONEP5*gamma/Sxx;
 628   B+=Db;
 629   A-=ONEP5*chi2/sum-xbar*Db;
 630   SB=sqrt((chi2/(n-2))/Sxx);
 631   SA=SB*sqrt(Sxx/sum+sqr(xbar));
 632   *a=A;
 633   *b=B;
 634   *sa=SA;
 635   *sb=SB;
 636 }
 637
 638