src/gromacs/statistics/statistics.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  * Copyright (c) 2012,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "statistics.h"
  40
  41 #include <math.h>
  42
  43 #include "gromacs/math/vec.h"
  44 #include "gromacs/utility/real.h"
  45 #include "gromacs/utility/smalloc.h"
  46
  47 static int gmx_dnint(double x)
  48 {
  49     return (int) (x+0.5);
  50 }
  51
  52 typedef struct gmx_stats {
  53     double  aa, a, b, sigma_aa, sigma_a, sigma_b, aver, sigma_aver, error;
  54     double  rmsd, Rdata, Rfit, Rfitaa, chi2, chi2aa;
  55     double *x, *y, *dx, *dy;
  56     int     computed;
  57     int     np, np_c, nalloc;
  58 } gmx_stats;
  59
  60 gmx_stats_t gmx_stats_init()
  61 {
  62     gmx_stats *stats;
  63
  64     snew(stats, 1);
  65
  66     return (gmx_stats_t) stats;
  67 }
  68
  69 int gmx_stats_get_npoints(gmx_stats_t gstats, int *N)
  70 {
  71     gmx_stats *stats = (gmx_stats *) gstats;
  72
  73     *N = stats->np;
  74
  75     return estatsOK;
  76 }
  77
  78 int gmx_stats_done(gmx_stats_t gstats)
  79 {
  80     gmx_stats *stats = (gmx_stats *) gstats;
  81
  82     sfree(stats->x);
  83     stats->x = NULL;
  84     sfree(stats->y);
  85     stats->y = NULL;
  86     sfree(stats->dx);
  87     stats->dx = NULL;
  88     sfree(stats->dy);
  89     stats->dy = NULL;
  90
  91     return estatsOK;
  92 }
  93
  94 int gmx_stats_add_point(gmx_stats_t gstats, double x, double y,
  95                         double dx, double dy)
  96 {
  97     gmx_stats *stats = (gmx_stats *) gstats;
  98     int        i;
  99
 100     if (stats->np+1 >= stats->nalloc)
 101     {
 102         if (stats->nalloc == 0)
 103         {
 104             stats->nalloc = 1024;
 105         }
 106         else
 107         {
 108             stats->nalloc *= 2;
 109         }
 110         srenew(stats->x, stats->nalloc);
 111         srenew(stats->y, stats->nalloc);
 112         srenew(stats->dx, stats->nalloc);
 113         srenew(stats->dy, stats->nalloc);
 114         for (i = stats->np; (i < stats->nalloc); i++)
 115         {
 116             stats->x[i]  = 0;
 117             stats->y[i]  = 0;
 118             stats->dx[i] = 0;
 119             stats->dy[i] = 0;
 120         }
 121     }
 122     stats->x[stats->np]  = x;
 123     stats->y[stats->np]  = y;
 124     stats->dx[stats->np] = dx;
 125     stats->dy[stats->np] = dy;
 126     stats->np++;
 127     stats->computed = 0;
 128
 129     return estatsOK;
 130 }
 131
 132 int gmx_stats_get_point(gmx_stats_t gstats, real *x, real *y,
 133                         real *dx, real *dy, real level)
 134 {
 135     gmx_stats *stats = (gmx_stats *) gstats;
 136     int        ok, outlier;
 137     real       rmsd, r;
 138
 139     if ((ok = gmx_stats_get_rmsd(gstats, &rmsd)) != estatsOK)
 140     {
 141         return ok;
 142     }
 143     outlier = 0;
 144     while ((outlier == 0) && (stats->np_c < stats->np))
 145     {
 146         r       = fabs(stats->x[stats->np_c] - stats->y[stats->np_c]);
 147         outlier = (r > rmsd*level);
 148         if (outlier)
 149         {
 150             if (NULL != x)
 151             {
 152                 *x  = stats->x[stats->np_c];
 153             }
 154             if (NULL != y)
 155             {
 156                 *y  = stats->y[stats->np_c];
 157             }
 158             if (NULL != dx)
 159             {
 160                 *dx = stats->dx[stats->np_c];
 161             }
 162             if (NULL != dy)
 163             {
 164                 *dy = stats->dy[stats->np_c];
 165             }
 166         }
 167         stats->np_c++;
 168
 169         if (outlier)
 170         {
 171             return estatsOK;
 172         }
 173     }
 174
 175     stats->np_c = 0;
 176
 177     return estatsNO_POINTS;
 178 }
 179
 180 int gmx_stats_add_points(gmx_stats_t gstats, int n, real *x, real *y,
 181                          real *dx, real *dy)
 182 {
 183     int i, ok;
 184
 185     for (i = 0; (i < n); i++)
 186     {
 187         if ((ok = gmx_stats_add_point(gstats, x[i], y[i],
 188                                       (NULL != dx) ? dx[i] : 0,
 189                                       (NULL != dy) ? dy[i] : 0)) != estatsOK)
 190         {
 191             return ok;
 192         }
 193     }
 194     return estatsOK;
 195 }
 196
 197 static int gmx_stats_compute(gmx_stats *stats, int weight)
 198 {
 199     double yy, yx, xx, sx, sy, dy, chi2, chi2aa, d2;
 200     double ssxx, ssyy, ssxy;
 201     double w, wtot, yx_nw, sy_nw, sx_nw, yy_nw, xx_nw, dx2, dy2;
 202     int    i, N;
 203
 204     N = stats->np;
 205     if (stats->computed == 0)
 206     {
 207         if (N < 1)
 208         {
 209             return estatsNO_POINTS;
 210         }
 211
 212         xx   = xx_nw = 0;
 213         yy   = yy_nw = 0;
 214         yx   = yx_nw = 0;
 215         sx   = sx_nw = 0;
 216         sy   = sy_nw = 0;
 217         wtot = 0;
 218         d2   = 0;
 219         for (i = 0; (i < N); i++)
 220         {
 221             d2 += dsqr(stats->x[i]-stats->y[i]);
 222             if ((stats->dy[i]) && (weight == elsqWEIGHT_Y))
 223             {
 224                 w = 1/dsqr(stats->dy[i]);
 225             }
 226             else
 227             {
 228                 w = 1;
 229             }
 230
 231             wtot  += w;
 232
 233             xx    += w*dsqr(stats->x[i]);
 234             xx_nw += dsqr(stats->x[i]);
 235
 236             yy    += w*dsqr(stats->y[i]);
 237             yy_nw += dsqr(stats->y[i]);
 238
 239             yx    += w*stats->y[i]*stats->x[i];
 240             yx_nw += stats->y[i]*stats->x[i];
 241
 242             sx    += w*stats->x[i];
 243             sx_nw += stats->x[i];
 244
 245             sy    += w*stats->y[i];
 246             sy_nw += stats->y[i];
 247         }
 248
 249         /* Compute average, sigma and error */
 250         stats->aver       = sy_nw/N;
 251         stats->sigma_aver = sqrt(yy_nw/N - dsqr(sy_nw/N));
 252         stats->error      = stats->sigma_aver/sqrt(N);
 253
 254         /* Compute RMSD between x and y */
 255         stats->rmsd = sqrt(d2/N);
 256
 257         /* Correlation coefficient for data */
 258         yx_nw       /= N;
 259         xx_nw       /= N;
 260         yy_nw       /= N;
 261         sx_nw       /= N;
 262         sy_nw       /= N;
 263         ssxx         = N*(xx_nw - dsqr(sx_nw));
 264         ssyy         = N*(yy_nw - dsqr(sy_nw));
 265         ssxy         = N*(yx_nw - (sx_nw*sy_nw));
 266         stats->Rdata = sqrt(dsqr(ssxy)/(ssxx*ssyy));
 267
 268         /* Compute straight line through datapoints, either with intercept
 269            zero (result in aa) or with intercept variable (results in a
 270            and b) */
 271         yx = yx/wtot;
 272         xx = xx/wtot;
 273         sx = sx/wtot;
 274         sy = sy/wtot;
 275
 276         stats->aa = (yx/xx);
 277         stats->a  = (yx-sx*sy)/(xx-sx*sx);
 278         stats->b  = (sy)-(stats->a)*(sx);
 279
 280         /* Compute chi2, deviation from a line y = ax+b. Also compute
 281            chi2aa which returns the deviation from a line y = ax. */
 282         chi2   = 0;
 283         chi2aa = 0;
 284         for (i = 0; (i < N); i++)
 285         {
 286             if (stats->dy[i] > 0)
 287             {
 288                 dy = stats->dy[i];
 289             }
 290             else
 291             {
 292                 dy = 1;
 293             }
 294             chi2aa += dsqr((stats->y[i]-(stats->aa*stats->x[i]))/dy);
 295             chi2   += dsqr((stats->y[i]-(stats->a*stats->x[i]+stats->b))/dy);
 296         }
 297         if (N > 2)
 298         {
 299             stats->chi2   = sqrt(chi2/(N-2));
 300             stats->chi2aa = sqrt(chi2aa/(N-2));
 301
 302             /* Look up equations! */
 303             dx2            = (xx-sx*sx);
 304             dy2            = (yy-sy*sy);
 305             stats->sigma_a = sqrt(stats->chi2/((N-2)*dx2));
 306             stats->sigma_b = stats->sigma_a*sqrt(xx);
 307             stats->Rfit    = fabs(ssxy)/sqrt(ssxx*ssyy);
 308             /*stats->a*sqrt(dx2/dy2);*/
 309             stats->Rfitaa  = stats->aa*sqrt(dx2/dy2);
 310         }
 311         else
 312         {
 313             stats->chi2    = 0;
 314             stats->chi2aa  = 0;
 315             stats->sigma_a = 0;
 316             stats->sigma_b = 0;
 317             stats->Rfit    = 0;
 318             stats->Rfitaa  = 0;
 319         }
 320
 321         stats->computed = 1;
 322     }
 323
 324     return estatsOK;
 325 }
 326
 327 int gmx_stats_get_ab(gmx_stats_t gstats, int weight,
 328                      real *a, real *b, real *da, real *db,
 329                      real *chi2, real *Rfit)
 330 {
 331     gmx_stats *stats = (gmx_stats *) gstats;
 332     int        ok;
 333
 334     if ((ok = gmx_stats_compute(stats, weight)) != estatsOK)
 335     {
 336         return ok;
 337     }
 338     if (NULL != a)
 339     {
 340         *a    = stats->a;
 341     }
 342     if (NULL != b)
 343     {
 344         *b    = stats->b;
 345     }
 346     if (NULL != da)
 347     {
 348         *da   = stats->sigma_a;
 349     }
 350     if (NULL != db)
 351     {
 352         *db   = stats->sigma_b;
 353     }
 354     if (NULL != chi2)
 355     {
 356         *chi2 = stats->chi2;
 357     }
 358     if (NULL != Rfit)
 359     {
 360         *Rfit = stats->Rfit;
 361     }
 362
 363     return estatsOK;
 364 }
 365
 366 int gmx_stats_get_a(gmx_stats_t gstats, int weight, real *a, real *da,
 367                     real *chi2, real *Rfit)
 368 {
 369     gmx_stats *stats = (gmx_stats *) gstats;
 370     int        ok;
 371
 372     if ((ok = gmx_stats_compute(stats, weight)) != estatsOK)
 373     {
 374         return ok;
 375     }
 376     if (NULL != a)
 377     {
 378         *a    = stats->aa;
 379     }
 380     if (NULL != da)
 381     {
 382         *da   = stats->sigma_aa;
 383     }
 384     if (NULL != chi2)
 385     {
 386         *chi2 = stats->chi2aa;
 387     }
 388     if (NULL != Rfit)
 389     {
 390         *Rfit = stats->Rfitaa;
 391     }
 392
 393     return estatsOK;
 394 }
 395
 396 int gmx_stats_get_average(gmx_stats_t gstats, real *aver)
 397 {
 398     gmx_stats *stats = (gmx_stats *) gstats;
 399     int        ok;
 400
 401     if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != estatsOK)
 402     {
 403         return ok;
 404     }
 405
 406     *aver = stats->aver;
 407
 408     return estatsOK;
 409 }
 410
 411 int gmx_stats_get_ase(gmx_stats_t gstats, real *aver, real *sigma, real *error)
 412 {
 413     gmx_stats *stats = (gmx_stats *) gstats;
 414     int        ok;
 415
 416     if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != estatsOK)
 417     {
 418         return ok;
 419     }
 420
 421     if (NULL != aver)
 422     {
 423         *aver  = stats->aver;
 424     }
 425     if (NULL != sigma)
 426     {
 427         *sigma = stats->sigma_aver;
 428     }
 429     if (NULL != error)
 430     {
 431         *error = stats->error;
 432     }
 433
 434     return estatsOK;
 435 }
 436
 437 int gmx_stats_get_sigma(gmx_stats_t gstats, real *sigma)
 438 {
 439     gmx_stats *stats = (gmx_stats *) gstats;
 440     int        ok;
 441
 442     if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != estatsOK)
 443     {
 444         return ok;
 445     }
 446
 447     *sigma = stats->sigma_aver;
 448
 449     return estatsOK;
 450 }
 451
 452 int gmx_stats_get_error(gmx_stats_t gstats, real *error)
 453 {
 454     gmx_stats *stats = (gmx_stats *) gstats;
 455     int        ok;
 456
 457     if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != estatsOK)
 458     {
 459         return ok;
 460     }
 461
 462     *error = stats->error;
 463
 464     return estatsOK;
 465 }
 466
 467 int gmx_stats_get_corr_coeff(gmx_stats_t gstats, real *R)
 468 {
 469     gmx_stats *stats = (gmx_stats *) gstats;
 470     int        ok;
 471
 472     if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != estatsOK)
 473     {
 474         return ok;
 475     }
 476
 477     *R = stats->Rdata;
 478
 479     return estatsOK;
 480 }
 481
 482 int gmx_stats_get_rmsd(gmx_stats_t gstats, real *rmsd)
 483 {
 484     gmx_stats *stats = (gmx_stats *) gstats;
 485     int        ok;
 486
 487     if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != estatsOK)
 488     {
 489         return ok;
 490     }
 491
 492     *rmsd = stats->rmsd;
 493
 494     return estatsOK;
 495 }
 496
 497 int gmx_stats_dump_xy(gmx_stats_t gstats, FILE *fp)
 498 {
 499     gmx_stats *stats = (gmx_stats *) gstats;
 500     int        i, ok;
 501
 502     for (i = 0; (i < stats->np); i++)
 503     {
 504         fprintf(fp, "%12g  %12g  %12g  %12g\n", stats->x[i], stats->y[i],
 505                 stats->dx[i], stats->dy[i]);
 506     }
 507
 508     return estatsOK;
 509 }
 510
 511 int gmx_stats_remove_outliers(gmx_stats_t gstats, double level)
 512 {
 513     gmx_stats *stats = (gmx_stats *) gstats;
 514     int        i, iter = 1, done = 0, ok;
 515     real       rmsd, r;
 516
 517     while ((stats->np >= 10) && !done)
 518     {
 519         if ((ok = gmx_stats_get_rmsd(gstats, &rmsd)) != estatsOK)
 520         {
 521             return ok;
 522         }
 523         done = 1;
 524         for (i = 0; (i < stats->np); )
 525         {
 526             r = fabs(stats->x[i]-stats->y[i]);
 527             if (r > level*rmsd)
 528             {
 529                 fprintf(stderr, "Removing outlier, iter = %d, rmsd = %g, x = %g, y = %g\n",
 530                         iter, rmsd, stats->x[i], stats->y[i]);
 531                 if (i < stats->np-1)
 532                 {
 533                     stats->x[i]  = stats->x[stats->np-1];
 534                     stats->y[i]  = stats->y[stats->np-1];
 535                     stats->dx[i] = stats->dx[stats->np-1];
 536                     stats->dy[i] = stats->dy[stats->np-1];
 537                 }
 538                 stats->np--;
 539                 done = 0;
 540             }
 541             else
 542             {
 543                 i++;
 544             }
 545         }
 546         iter++;
 547     }
 548
 549     return estatsOK;
 550 }
 551
 552 int gmx_stats_make_histogram(gmx_stats_t gstats, real binwidth, int *nb,
 553                              int ehisto, int normalized, real **x, real **y)
 554 {
 555     gmx_stats *stats = (gmx_stats *) gstats;
 556     int        i, ok, index = 0, nbins = *nb, *nindex;
 557     double     minx, maxx, maxy, miny, delta, dd, minh;
 558
 559     if (((binwidth <= 0) && (nbins <= 0)) ||
 560         ((binwidth > 0) && (nbins > 0)))
 561     {
 562         return estatsINVALID_INPUT;
 563     }
 564     if (stats->np <= 2)
 565     {
 566         return estatsNO_POINTS;
 567     }
 568     minx = maxx = stats->x[0];
 569     miny = maxy = stats->y[0];
 570     for (i = 1; (i < stats->np); i++)
 571     {
 572         miny = (stats->y[i] < miny) ? stats->y[i] : miny;
 573         maxy = (stats->y[i] > maxy) ? stats->y[i] : maxy;
 574         minx = (stats->x[i] < minx) ? stats->x[i] : minx;
 575         maxx = (stats->x[i] > maxx) ? stats->x[i] : maxx;
 576     }
 577     if (ehisto == ehistoX)
 578     {
 579         delta = maxx-minx;
 580         minh  = minx;
 581     }
 582     else if (ehisto == ehistoY)
 583     {
 584         delta = maxy-miny;
 585         minh  = miny;
 586     }
 587     else
 588     {
 589         return estatsINVALID_INPUT;
 590     }
 591
 592     if (binwidth == 0)
 593     {
 594         binwidth = (delta)/nbins;
 595     }
 596     else
 597     {
 598         nbins = gmx_dnint((delta)/binwidth + 0.5);
 599     }
 600     snew(*x, nbins);
 601     snew(nindex, nbins);
 602     for (i = 0; (i < nbins); i++)
 603     {
 604         (*x)[i] = minh + binwidth*(i+0.5);
 605     }
 606     if (normalized == 0)
 607     {
 608         dd = 1;
 609     }
 610     else
 611     {
 612         dd = 1.0/(binwidth*stats->np);
 613     }
 614
 615     snew(*y, nbins);
 616     for (i = 0; (i < stats->np); i++)
 617     {
 618         if (ehisto == ehistoY)
 619         {
 620             index = (stats->y[i]-miny)/binwidth;
 621         }
 622         else if (ehisto == ehistoX)
 623         {
 624             index = (stats->x[i]-minx)/binwidth;
 625         }
 626         if (index < 0)
 627         {
 628             index = 0;
 629         }
 630         if (index > nbins-1)
 631         {
 632             index = nbins-1;
 633         }
 634         (*y)[index] += dd;
 635         nindex[index]++;
 636     }
 637     if (*nb == 0)
 638     {
 639         *nb = nbins;
 640     }
 641     for (i = 0; (i < nbins); i++)
 642     {
 643         if (nindex[i] > 0)
 644         {
 645             (*y)[i] /= nindex[i];
 646         }
 647     }
 648
 649     sfree(nindex);
 650
 651     return estatsOK;
 652 }
 653
 654 static const char *stats_error[estatsNR] =
 655 {
 656     "All well in STATS land",
 657     "No points",
 658     "Not enough memory",
 659     "Invalid histogram input",
 660     "Unknown error",
 661     "Not implemented yet"
 662 };
 663
 664 const char *gmx_stats_message(int estats)
 665 {
 666     if ((estats >= 0) && (estats < estatsNR))
 667     {
 668         return stats_error[estats];
 669     }
 670     else
 671     {
 672         return stats_error[estatsERROR];
 673     }
 674 }
 675
 676 /* Old convenience functions, should be merged with the core
 677    statistics above. */
 678 int lsq_y_ax(int n, real x[], real y[], real *a)
 679 {
 680     gmx_stats_t lsq = gmx_stats_init();
 681     int         ok;
 682     real        da, chi2, Rfit;
 683
 684     gmx_stats_add_points(lsq, n, x, y, 0, 0);
 685     if ((ok = gmx_stats_get_a(lsq, elsqWEIGHT_NONE, a, &da, &chi2, &Rfit)) != estatsOK)
 686     {
 687         return ok;
 688     }
 689
 690     /*  int    i;
 691         double xx,yx;
 692
 693         yx=xx=0.0;
 694         for (i=0; i<n; i++) {
 695         yx+=y[i]*x[i];
 696         xx+=x[i]*x[i];
 697         }
 698      * a=yx/xx;
 699      */
 700     return estatsOK;
 701 }
 702
 703 static int low_lsq_y_ax_b(int n, real *xr, double *xd, real yr[],
 704                           real *a, real *b, real *r, real *chi2)
 705 {
 706     int         i, ok;
 707     gmx_stats_t lsq;
 708
 709     lsq = gmx_stats_init();
 710     for (i = 0; (i < n); i++)
 711     {
 712         if ((ok = gmx_stats_add_point(lsq, (NULL != xd) ? xd[i] : xr[i], yr[i], 0, 0))
 713             != estatsOK)
 714         {
 715             return ok;
 716         }
 717     }
 718     if ((ok = gmx_stats_get_ab(lsq, elsqWEIGHT_NONE, a, b, NULL, NULL, chi2, r)) != estatsOK)
 719     {
 720         return ok;
 721     }
 722
 723     return estatsOK;
 724     /*
 725        double x,y,yx,xx,yy,sx,sy,chi2;
 726
 727        yx=xx=yy=sx=sy=0.0;
 728        for (i=0; i<n; i++) {
 729        if (xd != NULL) {
 730        x = xd[i];
 731        } else {
 732        x = xr[i];
 733        }
 734        y =   yr[i];
 735
 736        yx += y*x;
 737        xx += x*x;
 738        yy += y*y;
 739        sx += x;
 740        sy += y;
 741        }
 742      * a = (n*yx-sy*sx)/(n*xx-sx*sx);
 743      * b = (sy-(*a)*sx)/n;
 744      * r = sqrt((xx-sx*sx)/(yy-sy*sy));
 745
 746        chi2 = 0;
 747        if (xd != NULL) {
 748        for(i=0; i<n; i++)
 749        chi2 += dsqr(yr[i] - ((*a)*xd[i] + (*b)));
 750        } else {
 751        for(i=0; i<n; i++)
 752        chi2 += dsqr(yr[i] - ((*a)*xr[i] + (*b)));
 753        }
 754
 755        if (n > 2)
 756        return sqrt(chi2/(n-2));
 757        else
 758        return 0;
 759      */
 760 }
 761
 762 int lsq_y_ax_b(int n, real x[], real y[], real *a, real *b, real *r, real *chi2)
 763 {
 764     return low_lsq_y_ax_b(n, x, NULL, y, a, b, r, chi2);
 765 }
 766
 767 int lsq_y_ax_b_xdouble(int n, double x[], real y[], real *a, real *b,
 768                        real *r, real *chi2)
 769 {
 770     return low_lsq_y_ax_b(n, NULL, x, y, a, b, r, chi2);
 771 }
 772
 773 int lsq_y_ax_b_error(int n, real x[], real y[], real dy[],
 774                      real *a, real *b, real *da, real *db,
 775                      real *r, real *chi2)
 776 {
 777     gmx_stats_t lsq;
 778     int         i, ok;
 779
 780     lsq = gmx_stats_init();
 781     for (i = 0; (i < n); i++)
 782     {
 783         if ((ok = gmx_stats_add_point(lsq, x[i], y[i], 0, dy[i])) != estatsOK)
 784         {
 785             return ok;
 786         }
 787     }
 788     if ((ok = gmx_stats_get_ab(lsq, elsqWEIGHT_Y, a, b, da, db, chi2, r)) != estatsOK)
 789     {
 790         return ok;
 791     }
 792     if ((ok = gmx_stats_done(lsq)) != estatsOK)
 793     {
 794         return ok;
 795     }
 796     sfree(lsq);
 797
 798     return estatsOK;
 799     /*
 800        double sxy,sxx,syy,sx,sy,w,s_2,dx2,dy2,mins;
 801
 802        sxy=sxx=syy=sx=sy=w=0.0;
 803        mins = dy[0];
 804        for(i=1; (i<n); i++)
 805        mins = min(mins,dy[i]);
 806        if (mins <= 0)
 807        gmx_fatal(FARGS,"Zero or negative weigths in linear regression analysis");
 808
 809        for (i=0; i<n; i++) {
 810        s_2  = dsqr(1.0/dy[i]);
 811        sxx += s_2*dsqr(x[i]);
 812        sxy += s_2*y[i]*x[i];
 813        syy += s_2*dsqr(y[i]);
 814        sx  += s_2*x[i];
 815        sy  += s_2*y[i];
 816        w   += s_2;
 817        }
 818        sxx = sxx/w;
 819        sxy = sxy/w;
 820        syy = syy/w;
 821        sx  = sx/w;
 822        sy  = sy/w;
 823        dx2 = (sxx-sx*sx);
 824        dy2 = (syy-sy*sy);
 825      * a=(sxy-sy*sx)/dx2;
 826      * b=(sy-(*a)*sx);
 827
 828      * chi2=0;
 829        for(i=0; i<n; i++)
 830      * chi2+=dsqr((y[i]-((*a)*x[i]+(*b)))/dy[i]);
 831      * chi2 = *chi2/w;
 832
 833      * da = sqrt(*chi2/((n-2)*dx2));
 834      * db = *da*sqrt(sxx);
 835      * r  = *a*sqrt(dx2/dy2);
 836
 837        if (debug)
 838        fprintf(debug,"sx = %g, sy = %g, sxy = %g, sxx = %g, w = %g\n"
 839        "chi2 = %g, dx2 = %g\n",
 840        sx,sy,sxy,sxx,w,*chi2,dx2);
 841
 842        if (n > 2)
 843      * chi2 = sqrt(*chi2/(n-2));
 844        else
 845      * chi2 = 0;
 846      */
 847 }