htswanalysis/MACS/lib/gsl/gsl-1.11/multifit/multilinear.c

   1 /* multifit/multilinear.c
   2  *
   3  * Copyright (C) 2000, 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <gsl/gsl_errno.h>
  22 #include <gsl/gsl_multifit.h>
  23 #include <gsl/gsl_blas.h>
  24 #include <gsl/gsl_linalg.h>
  25
  26 /* Fit
  27  *
  28  *  y = X c
  29  *
  30  *  where X is an M x N matrix of M observations for N variables.
  31  *
  32  */
  33
  34 int
  35 gsl_multifit_linear (const gsl_matrix * X,
  36                      const gsl_vector * y,
  37                      gsl_vector * c,
  38                      gsl_matrix * cov,
  39                      double *chisq, gsl_multifit_linear_workspace * work)
  40 {
  41   size_t rank;
  42   int status  = gsl_multifit_linear_svd (X, y, GSL_DBL_EPSILON, &rank, c,
  43                                          cov, chisq, work);
  44   return status;
  45 }
  46
  47 /* Handle the general case of the SVD with tolerance and rank */
  48
  49 int
  50 gsl_multifit_linear_svd (const gsl_matrix * X,
  51                          const gsl_vector * y,
  52                          double tol,
  53                          size_t * rank,
  54                          gsl_vector * c,
  55                          gsl_matrix * cov,
  56                          double *chisq, gsl_multifit_linear_workspace * work)
  57 {
  58   if (X->size1 != y->size)
  59     {
  60       GSL_ERROR
  61         ("number of observations in y does not match rows of matrix X",
  62          GSL_EBADLEN);
  63     }
  64   else if (X->size2 != c->size)
  65     {
  66       GSL_ERROR ("number of parameters c does not match columns of matrix X",
  67                  GSL_EBADLEN);
  68     }
  69   else if (cov->size1 != cov->size2)
  70     {
  71       GSL_ERROR ("covariance matrix is not square", GSL_ENOTSQR);
  72     }
  73   else if (c->size != cov->size1)
  74     {
  75       GSL_ERROR
  76         ("number of parameters does not match size of covariance matrix",
  77          GSL_EBADLEN);
  78     }
  79   else if (X->size1 != work->n || X->size2 != work->p)
  80     {
  81       GSL_ERROR
  82         ("size of workspace does not match size of observation matrix",
  83          GSL_EBADLEN);
  84     }
  85   else if (tol <= 0)
  86     {
  87       GSL_ERROR ("tolerance must be positive", GSL_EINVAL);
  88     }
  89   else
  90     {
  91       const size_t n = X->size1;
  92       const size_t p = X->size2;
  93
  94       size_t i, j, p_eff;
  95
  96       gsl_matrix *A = work->A;
  97       gsl_matrix *Q = work->Q;
  98       gsl_matrix *QSI = work->QSI;
  99       gsl_vector *S = work->S;
 100       gsl_vector *xt = work->xt;
 101       gsl_vector *D = work->D;
 102
 103       /* Copy X to workspace,  A <= X */
 104
 105       gsl_matrix_memcpy (A, X);
 106
 107       /* Balance the columns of the matrix A */
 108
 109       gsl_linalg_balance_columns (A, D);
 110
 111       /* Decompose A into U S Q^T */
 112
 113       gsl_linalg_SV_decomp_mod (A, QSI, Q, S, xt);
 114
 115       /* Solve y = A c for c */
 116
 117       gsl_blas_dgemv (CblasTrans, 1.0, A, y, 0.0, xt);
 118
 119       /* Scale the matrix Q,  Q' = Q S^-1 */
 120
 121       gsl_matrix_memcpy (QSI, Q);
 122
 123       {
 124         double alpha0 = gsl_vector_get (S, 0);
 125         p_eff = 0;
 126
 127         for (j = 0; j < p; j++)
 128           {
 129             gsl_vector_view column = gsl_matrix_column (QSI, j);
 130             double alpha = gsl_vector_get (S, j);
 131
 132             if (alpha <= tol * alpha0) {
 133               alpha = 0.0;
 134             } else {
 135               alpha = 1.0 / alpha;
 136               p_eff++;
 137             }
 138
 139             gsl_vector_scale (&column.vector, alpha);
 140           }
 141
 142         *rank = p_eff;
 143       }
 144
 145       gsl_vector_set_zero (c);
 146
 147       gsl_blas_dgemv (CblasNoTrans, 1.0, QSI, xt, 0.0, c);
 148
 149       /* Unscale the balancing factors */
 150
 151       gsl_vector_div (c, D);
 152
 153       /* Compute chisq, from residual r = y - X c */
 154
 155       {
 156         double s2 = 0, r2 = 0;
 157
 158         for (i = 0; i < n; i++)
 159           {
 160             double yi = gsl_vector_get (y, i);
 161             gsl_vector_const_view row = gsl_matrix_const_row (X, i);
 162             double y_est, ri;
 163             gsl_blas_ddot (&row.vector, c, &y_est);
 164             ri = yi - y_est;
 165             r2 += ri * ri;
 166           }
 167
 168         s2 = r2 / (n - p_eff);   /* p_eff == rank */
 169
 170         *chisq = r2;
 171
 172         /* Form variance-covariance matrix cov = s2 * (Q S^-1) (Q S^-1)^T */
 173
 174         for (i = 0; i < p; i++)
 175           {
 176             gsl_vector_view row_i = gsl_matrix_row (QSI, i);
 177             double d_i = gsl_vector_get (D, i);
 178
 179             for (j = i; j < p; j++)
 180               {
 181                 gsl_vector_view row_j = gsl_matrix_row (QSI, j);
 182                 double d_j = gsl_vector_get (D, j);
 183                 double s;
 184
 185                 gsl_blas_ddot (&row_i.vector, &row_j.vector, &s);
 186
 187                 gsl_matrix_set (cov, i, j, s * s2 / (d_i * d_j));
 188                 gsl_matrix_set (cov, j, i, s * s2 / (d_i * d_j));
 189               }
 190           }
 191       }
 192
 193       return GSL_SUCCESS;
 194     }
 195 }
 196
 197 int
 198 gsl_multifit_wlinear (const gsl_matrix * X,
 199                       const gsl_vector * w,
 200                       const gsl_vector * y,
 201                       gsl_vector * c,
 202                       gsl_matrix * cov,
 203                       double *chisq, gsl_multifit_linear_workspace * work)
 204 {
 205   size_t rank;
 206   int status  = gsl_multifit_wlinear_svd (X, w, y, GSL_DBL_EPSILON, &rank, c,
 207                                           cov, chisq, work);
 208   return status;
 209 }
 210
 211
 212 int
 213 gsl_multifit_wlinear_svd (const gsl_matrix * X,
 214                           const gsl_vector * w,
 215                           const gsl_vector * y,
 216                           double tol,
 217                           size_t * rank,
 218                           gsl_vector * c,
 219                           gsl_matrix * cov,
 220                           double *chisq, gsl_multifit_linear_workspace * work)
 221 {
 222   if (X->size1 != y->size)
 223     {
 224       GSL_ERROR
 225         ("number of observations in y does not match rows of matrix X",
 226          GSL_EBADLEN);
 227     }
 228   else if (X->size2 != c->size)
 229     {
 230       GSL_ERROR ("number of parameters c does not match columns of matrix X",
 231                  GSL_EBADLEN);
 232     }
 233   else if (w->size != y->size)
 234     {
 235       GSL_ERROR ("number of weights does not match number of observations",
 236                  GSL_EBADLEN);
 237     }
 238   else if (cov->size1 != cov->size2)
 239     {
 240       GSL_ERROR ("covariance matrix is not square", GSL_ENOTSQR);
 241     }
 242   else if (c->size != cov->size1)
 243     {
 244       GSL_ERROR
 245         ("number of parameters does not match size of covariance matrix",
 246          GSL_EBADLEN);
 247     }
 248   else if (X->size1 != work->n || X->size2 != work->p)
 249     {
 250       GSL_ERROR
 251         ("size of workspace does not match size of observation matrix",
 252          GSL_EBADLEN);
 253     }
 254   else
 255     {
 256       const size_t n = X->size1;
 257       const size_t p = X->size2;
 258
 259       size_t i, j, p_eff;
 260
 261       gsl_matrix *A = work->A;
 262       gsl_matrix *Q = work->Q;
 263       gsl_matrix *QSI = work->QSI;
 264       gsl_vector *S = work->S;
 265       gsl_vector *t = work->t;
 266       gsl_vector *xt = work->xt;
 267       gsl_vector *D = work->D;
 268
 269       /* Scale X,  A = sqrt(w) X */
 270
 271       gsl_matrix_memcpy (A, X);
 272
 273       for (i = 0; i < n; i++)
 274         {
 275           double wi = gsl_vector_get (w, i);
 276
 277           if (wi < 0)
 278             wi = 0;
 279
 280           {
 281             gsl_vector_view row = gsl_matrix_row (A, i);
 282             gsl_vector_scale (&row.vector, sqrt (wi));
 283           }
 284         }
 285
 286       /* Balance the columns of the matrix A */
 287
 288       gsl_linalg_balance_columns (A, D);
 289
 290       /* Decompose A into U S Q^T */
 291
 292       gsl_linalg_SV_decomp_mod (A, QSI, Q, S, xt);
 293
 294       /* Solve sqrt(w) y = A c for c, by first computing t = sqrt(w) y */
 295
 296       for (i = 0; i < n; i++)
 297         {
 298           double wi = gsl_vector_get (w, i);
 299           double yi = gsl_vector_get (y, i);
 300           if (wi < 0)
 301             wi = 0;
 302           gsl_vector_set (t, i, sqrt (wi) * yi);
 303         }
 304
 305       gsl_blas_dgemv (CblasTrans, 1.0, A, t, 0.0, xt);
 306
 307       /* Scale the matrix Q,  Q' = Q S^-1 */
 308
 309       gsl_matrix_memcpy (QSI, Q);
 310
 311       {
 312         double alpha0 = gsl_vector_get (S, 0);
 313         p_eff = 0;
 314
 315         for (j = 0; j < p; j++)
 316           {
 317             gsl_vector_view column = gsl_matrix_column (QSI, j);
 318             double alpha = gsl_vector_get (S, j);
 319
 320             if (alpha <= tol * alpha0) {
 321               alpha = 0.0;
 322             } else {
 323               alpha = 1.0 / alpha;
 324               p_eff++;
 325             }
 326
 327             gsl_vector_scale (&column.vector, alpha);
 328           }
 329
 330         *rank = p_eff;
 331       }
 332
 333       gsl_vector_set_zero (c);
 334
 335       /* Solution */
 336
 337       gsl_blas_dgemv (CblasNoTrans, 1.0, QSI, xt, 0.0, c);
 338
 339       /* Unscale the balancing factors */
 340
 341       gsl_vector_div (c, D);
 342
 343       /* Form covariance matrix cov = (Q S^-1) (Q S^-1)^T */
 344
 345       for (i = 0; i < p; i++)
 346         {
 347           gsl_vector_view row_i = gsl_matrix_row (QSI, i);
 348           double d_i = gsl_vector_get (D, i);
 349
 350           for (j = i; j < p; j++)
 351             {
 352               gsl_vector_view row_j = gsl_matrix_row (QSI, j);
 353               double d_j = gsl_vector_get (D, j);
 354               double s;
 355
 356               gsl_blas_ddot (&row_i.vector, &row_j.vector, &s);
 357
 358               gsl_matrix_set (cov, i, j, s / (d_i * d_j));
 359               gsl_matrix_set (cov, j, i, s / (d_i * d_j));
 360             }
 361         }
 362
 363       /* Compute chisq, from residual r = y - X c */
 364
 365       {
 366         double r2 = 0;
 367
 368         for (i = 0; i < n; i++)
 369           {
 370             double yi = gsl_vector_get (y, i);
 371             double wi = gsl_vector_get (w, i);
 372             gsl_vector_const_view row = gsl_matrix_const_row (X, i);
 373             double y_est, ri;
 374             gsl_blas_ddot (&row.vector, c, &y_est);
 375             ri = yi - y_est;
 376             r2 += wi * ri * ri;
 377           }
 378
 379         *chisq = r2;
 380       }
 381
 382       return GSL_SUCCESS;
 383     }
 384 }
 385
 386
 387 int
 388 gsl_multifit_linear_est (const gsl_vector * x,
 389                          const gsl_vector * c,
 390                          const gsl_matrix * cov, double *y, double *y_err)
 391 {
 392
 393   if (x->size != c->size)
 394     {
 395       GSL_ERROR ("number of parameters c does not match number of observations x",
 396          GSL_EBADLEN);
 397     }
 398   else if (cov->size1 != cov->size2)
 399     {
 400       GSL_ERROR ("covariance matrix is not square", GSL_ENOTSQR);
 401     }
 402   else if (c->size != cov->size1)
 403     {
 404       GSL_ERROR ("number of parameters c does not match size of covariance matrix cov",
 405          GSL_EBADLEN);
 406     }
 407   else
 408     {
 409       size_t i, j;
 410       double var = 0;
 411
 412       gsl_blas_ddot(x, c, y);       /* y = x.c */
 413
 414       /* var = x' cov x */
 415
 416       for (i = 0; i < x->size; i++)
 417         {
 418           const double xi = gsl_vector_get (x, i);
 419           var += xi * xi * gsl_matrix_get (cov, i, i);
 420
 421           for (j = 0; j < i; j++)
 422             {
 423               const double xj = gsl_vector_get (x, j);
 424               var += 2 * xi * xj * gsl_matrix_get (cov, i, j);
 425             }
 426         }
 427
 428       *y_err = sqrt (var);
 429
 430       return GSL_SUCCESS;
 431     }
 432 }
 433
 434 /*
 435 gsl_multifit_linear_residuals()
 436   Compute vector of residuals from fit
 437
 438 Inputs: X - design matrix
 439         y - rhs vector
 440         c - fit coefficients
 441         r - (output) where to store residuals
 442 */
 443
 444 int
 445 gsl_multifit_linear_residuals (const gsl_matrix *X, const gsl_vector *y,
 446                                const gsl_vector *c, gsl_vector *r)
 447 {
 448   if (X->size1 != y->size)
 449     {
 450       GSL_ERROR
 451         ("number of observations in y does not match rows of matrix X",
 452          GSL_EBADLEN);
 453     }
 454   else if (X->size2 != c->size)
 455     {
 456       GSL_ERROR ("number of parameters c does not match columns of matrix X",
 457                  GSL_EBADLEN);
 458     }
 459   else if (y->size != r->size)
 460     {
 461       GSL_ERROR ("number of observations in y does not match number of residuals",
 462                  GSL_EBADLEN);
 463     }
 464   else
 465     {
 466       size_t i;
 467
 468       for (i = 0; i < y->size; ++i)
 469         {
 470           double yi = gsl_vector_get(y, i);
 471           gsl_vector_const_view row = gsl_matrix_const_row(X, i);
 472           double y_est, ri;
 473
 474           gsl_blas_ddot(&row.vector, c, &y_est);
 475           ri = yi - y_est;
 476
 477           gsl_vector_set(r, i, ri);
 478         }
 479
 480       return GSL_SUCCESS;
 481     }
 482 } /* gsl_multifit_linear_residuals() */