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