include/gsl/gsl_multifit.h

0001 /* multifit/gsl_multifit.h
0002  *
0003  * Copyright (C) 2000, 2007, 2010 Brian Gough
0004  * Copyright (C) 2013, Patrick Alken
0005  *
0006  * This program is free software; you can redistribute it and/or modify
0007  * it under the terms of the GNU General Public License as published by
0008  * the Free Software Foundation; either version 3 of the License, or (at
0009  * your option) any later version.
0010  *
0011  * This program is distributed in the hope that it will be useful, but
0012  * WITHOUT ANY WARRANTY; without even the implied warranty of
0013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0014  * General Public License for more details.
0015  *
0016  * You should have received a copy of the GNU General Public License
0017  * along with this program; if not, write to the Free Software
0018  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
0019  */
0020
0021 #ifndef __GSL_MULTIFIT_H__
0022 #define __GSL_MULTIFIT_H__
0023
0024 #include <stdlib.h>
0025 #include <gsl/gsl_math.h>
0026 #include <gsl/gsl_vector.h>
0027 #include <gsl/gsl_matrix.h>
0028 #include <gsl/gsl_types.h>
0029
0030 #undef __BEGIN_DECLS
0031 #undef __END_DECLS
0032 #ifdef __cplusplus
0033 # define __BEGIN_DECLS extern "C" {
0034 # define __END_DECLS }
0035 #else
0036 # define __BEGIN_DECLS /* empty */
0037 # define __END_DECLS /* empty */
0038 #endif
0039
0040 __BEGIN_DECLS
0041
0042 typedef struct
0043 {
0044   size_t nmax;         /* maximum number of observations */
0045   size_t pmax;         /* maximum number of parameters */
0046   size_t n;            /* number of observations in current SVD decomposition */
0047   size_t p;            /* number of parameters in current SVD decomposition */
0048   gsl_matrix * A;      /* least squares matrix for SVD, n-by-p */
0049   gsl_matrix * Q;
0050   gsl_matrix * QSI;
0051   gsl_vector * S;
0052   gsl_vector * t;
0053   gsl_vector * xt;
0054   gsl_vector * D;
0055   double rcond;        /* reciprocal condition number */
0056 }
0057 gsl_multifit_linear_workspace;
0058
0059 gsl_multifit_linear_workspace *
0060 gsl_multifit_linear_alloc (const size_t n, const size_t p);
0061
0062 void
0063 gsl_multifit_linear_free (gsl_multifit_linear_workspace * w);
0064
0065 int
0066 gsl_multifit_linear (const gsl_matrix * X,
0067                      const gsl_vector * y,
0068                      gsl_vector * c,
0069                      gsl_matrix * cov,
0070                      double * chisq,
0071                      gsl_multifit_linear_workspace * work);
0072
0073 int
0074 gsl_multifit_linear_tsvd (const gsl_matrix * X,
0075                           const gsl_vector * y,
0076                           const double tol,
0077                           gsl_vector * c,
0078                           gsl_matrix * cov,
0079                           double * chisq,
0080                           size_t * rank,
0081                           gsl_multifit_linear_workspace * work);
0082
0083 int
0084 gsl_multifit_linear_svd (const gsl_matrix * X,
0085                          gsl_multifit_linear_workspace * work);
0086
0087 int
0088 gsl_multifit_linear_bsvd (const gsl_matrix * X,
0089                           gsl_multifit_linear_workspace * work);
0090
0091 size_t
0092 gsl_multifit_linear_rank(const double tol, const gsl_multifit_linear_workspace * work);
0093
0094 int
0095 gsl_multifit_linear_solve (const double lambda,
0096                            const gsl_matrix * X,
0097                            const gsl_vector * y,
0098                            gsl_vector * c,
0099                            double *rnorm,
0100                            double *snorm,
0101                            gsl_multifit_linear_workspace * work);
0102
0103 int
0104 gsl_multifit_linear_applyW(const gsl_matrix * X,
0105                            const gsl_vector * w,
0106                            const gsl_vector * y,
0107                            gsl_matrix * WX,
0108                            gsl_vector * Wy);
0109
0110 int
0111 gsl_multifit_linear_stdform1 (const gsl_vector * L,
0112                               const gsl_matrix * X,
0113                               const gsl_vector * y,
0114                               gsl_matrix * Xs,
0115                               gsl_vector * ys,
0116                               gsl_multifit_linear_workspace * work);
0117
0118 int
0119 gsl_multifit_linear_wstdform1 (const gsl_vector * L,
0120                                const gsl_matrix * X,
0121                                const gsl_vector * w,
0122                                const gsl_vector * y,
0123                                gsl_matrix * Xs,
0124                                gsl_vector * ys,
0125                                gsl_multifit_linear_workspace * work);
0126
0127 int
0128 gsl_multifit_linear_L_decomp (gsl_matrix * L, gsl_vector * tau);
0129
0130 int
0131 gsl_multifit_linear_stdform2 (const gsl_matrix * LQR,
0132                               const gsl_vector * Ltau,
0133                               const gsl_matrix * X,
0134                               const gsl_vector * y,
0135                               gsl_matrix * Xs,
0136                               gsl_vector * ys,
0137                               gsl_matrix * M,
0138                               gsl_multifit_linear_workspace * work);
0139
0140 int
0141 gsl_multifit_linear_wstdform2 (const gsl_matrix * LQR,
0142                                const gsl_vector * Ltau,
0143                                const gsl_matrix * X,
0144                                const gsl_vector * w,
0145                                const gsl_vector * y,
0146                                gsl_matrix * Xs,
0147                                gsl_vector * ys,
0148                                gsl_matrix * M,
0149                                gsl_multifit_linear_workspace * work);
0150
0151 int
0152 gsl_multifit_linear_genform1 (const gsl_vector * L,
0153                               const gsl_vector * cs,
0154                               gsl_vector * c,
0155                               gsl_multifit_linear_workspace * work);
0156
0157 int
0158 gsl_multifit_linear_genform2 (const gsl_matrix * LQR,
0159                               const gsl_vector * Ltau,
0160                               const gsl_matrix * X,
0161                               const gsl_vector * y,
0162                               const gsl_vector * cs,
0163                               const gsl_matrix * M,
0164                               gsl_vector * c,
0165                               gsl_multifit_linear_workspace * work);
0166
0167 int
0168 gsl_multifit_linear_wgenform2 (const gsl_matrix * LQR,
0169                                const gsl_vector * Ltau,
0170                                const gsl_matrix * X,
0171                                const gsl_vector * w,
0172                                const gsl_vector * y,
0173                                const gsl_vector * cs,
0174                                const gsl_matrix * M,
0175                                gsl_vector * c,
0176                                gsl_multifit_linear_workspace * work);
0177
0178 int
0179 gsl_multifit_linear_lreg (const double smin, const double smax,
0180                           gsl_vector * reg_param);
0181
0182 int
0183 gsl_multifit_linear_lcurve (const gsl_vector * y,
0184                             gsl_vector * reg_param,
0185                             gsl_vector * rho, gsl_vector * eta,
0186                             gsl_multifit_linear_workspace * work);
0187
0188 int
0189 gsl_multifit_linear_lcurvature (const gsl_vector * y,
0190                                 const gsl_vector * reg_param,
0191                                 const gsl_vector * rho,
0192                                 const gsl_vector * eta,
0193                                 gsl_vector * kappa,
0194                                 gsl_multifit_linear_workspace * work);
0195
0196 int
0197 gsl_multifit_linear_lcurvature_menger(const gsl_vector *rho,
0198                                       const gsl_vector *eta,
0199                                       gsl_vector *kappa);
0200
0201 int
0202 gsl_multifit_linear_lcorner(const gsl_vector *rho,
0203                             const gsl_vector *eta,
0204                             size_t *idx);
0205
0206 int
0207 gsl_multifit_linear_lcorner2(const gsl_vector *reg_param,
0208                              const gsl_vector *eta,
0209                              size_t *idx);
0210
0211 int
0212 gsl_multifit_linear_Lk(const size_t p, const size_t k, gsl_matrix *L);
0213
0214 int
0215 gsl_multifit_linear_Lsobolev(const size_t p, const size_t kmax,
0216                              const gsl_vector *alpha, gsl_matrix *L,
0217                              gsl_multifit_linear_workspace *work);
0218
0219 int
0220 gsl_multifit_wlinear (const gsl_matrix * X,
0221                       const gsl_vector * w,
0222                       const gsl_vector * y,
0223                       gsl_vector * c,
0224                       gsl_matrix * cov,
0225                       double * chisq,
0226                       gsl_multifit_linear_workspace * work);
0227
0228 int
0229 gsl_multifit_wlinear_tsvd (const gsl_matrix * X,
0230                            const gsl_vector * w,
0231                            const gsl_vector * y,
0232                            const double tol,
0233                            gsl_vector * c,
0234                            gsl_matrix * cov,
0235                            double * chisq,
0236                            size_t * rank,
0237                            gsl_multifit_linear_workspace * work);
0238
0239 int
0240 gsl_multifit_wlinear_svd (const gsl_matrix * X,
0241                           const gsl_vector * w,
0242                           const gsl_vector * y,
0243                           double tol,
0244                           size_t * rank,
0245                           gsl_vector * c,
0246                           gsl_matrix * cov,
0247                           double *chisq,
0248                           gsl_multifit_linear_workspace * work);
0249
0250 int
0251 gsl_multifit_wlinear_usvd (const gsl_matrix * X,
0252                            const gsl_vector * w,
0253                            const gsl_vector * y,
0254                            double tol,
0255                            size_t * rank,
0256                            gsl_vector * c,
0257                            gsl_matrix * cov,
0258                            double *chisq,
0259                            gsl_multifit_linear_workspace * work);
0260
0261 int
0262 gsl_multifit_linear_est (const gsl_vector * x,
0263                          const gsl_vector * c,
0264                          const gsl_matrix * cov, double *y, double *y_err);
0265
0266 double
0267 gsl_multifit_linear_rcond (const gsl_multifit_linear_workspace * w);
0268
0269 int
0270 gsl_multifit_linear_residuals (const gsl_matrix *X, const gsl_vector *y,
0271                                const gsl_vector *c, gsl_vector *r);
0272
0273 /* gcv.c */
0274 int
0275 gsl_multifit_linear_gcv_init(const gsl_vector * y,
0276                              gsl_vector * reg_param,
0277                              gsl_vector * UTy,
0278                              double * delta0,
0279                              gsl_multifit_linear_workspace * work);
0280
0281 int
0282 gsl_multifit_linear_gcv_curve(const gsl_vector * reg_param,
0283                               const gsl_vector * UTy,
0284                               const double delta0,
0285                               gsl_vector * G,
0286                               gsl_multifit_linear_workspace * work);
0287
0288 int
0289 gsl_multifit_linear_gcv_min(const gsl_vector * reg_param,
0290                             const gsl_vector * UTy,
0291                             const gsl_vector * G,
0292                             const double delta0,
0293                             double * lambda,
0294                             gsl_multifit_linear_workspace * work);
0295
0296 double
0297 gsl_multifit_linear_gcv_calc(const double lambda,
0298                              const gsl_vector * UTy,
0299                              const double delta0,
0300                              gsl_multifit_linear_workspace * work);
0301
0302 int
0303 gsl_multifit_linear_gcv(const gsl_vector * y,
0304                         gsl_vector * reg_param,
0305                         gsl_vector * G,
0306                         double * lambda,
0307                         double * G_lambda,
0308                         gsl_multifit_linear_workspace * work);
0309
0310 typedef struct
0311 {
0312   const char * name;     /* method name */
0313   int (*wfun)(const gsl_vector *r, gsl_vector *w);
0314   int (*psi_deriv)(const gsl_vector *r, gsl_vector *dpsi);
0315   double tuning_default; /* default tuning constant */
0316 } gsl_multifit_robust_type;
0317
0318 typedef struct
0319 {
0320   double sigma_ols;    /* OLS estimate of sigma */
0321   double sigma_mad;    /* MAD estimate of sigma */
0322   double sigma_rob;    /* robust estimate of sigma */
0323   double sigma;        /* final estimate of sigma */
0324   double Rsq;          /* R^2 coefficient of determination */
0325   double adj_Rsq;      /* degree of freedom adjusted R^2 */
0326   double rmse;         /* root mean squared error */
0327   double sse;          /* residual sum of squares */
0328   size_t dof;          /* degrees of freedom */
0329   size_t numit;        /* number of iterations */
0330   gsl_vector *weights; /* final weights */
0331   gsl_vector *r;       /* final residuals y - X c */
0332 } gsl_multifit_robust_stats;
0333
0334 typedef struct
0335 {
0336   size_t n;            /* number of observations */
0337   size_t p;            /* number of parameters */
0338   size_t numit;        /* number of iterations */
0339   size_t maxiter;      /* maximum iterations */
0340   const gsl_multifit_robust_type *type;
0341   double tune;         /* tuning parameter */
0342
0343   gsl_vector *r;       /* residuals at current iteration */
0344   gsl_vector *weights; /* weights at current iteration */
0345   gsl_vector *c_prev;  /* coefficients from previous iteration */
0346   gsl_vector *resfac;  /* multiplicative factors for residuals */
0347
0348   gsl_vector *psi;     /* psi(r) */
0349   gsl_vector *dpsi;    /* psi'(r) */
0350
0351   gsl_matrix *QSI;     /* Q S^{-1} of original matrix X */
0352   gsl_vector *D;       /* balancing parameters of original matrix X */
0353
0354   gsl_vector *workn;   /* workspace of length n */
0355
0356   gsl_multifit_robust_stats stats; /* various statistics */
0357
0358   gsl_multifit_linear_workspace *multifit_p;
0359 } gsl_multifit_robust_workspace;
0360
0361 /* available types */
0362 GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_default;
0363 GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_bisquare;
0364 GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_cauchy;
0365 GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_fair;
0366 GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_huber;
0367 GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_ols;
0368 GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_welsch;
0369
0370 gsl_multifit_robust_workspace *gsl_multifit_robust_alloc(const gsl_multifit_robust_type *T,
0371                                                          const size_t n, const size_t p);
0372 void gsl_multifit_robust_free(gsl_multifit_robust_workspace *w);
0373 int gsl_multifit_robust_tune(const double tune,
0374                              gsl_multifit_robust_workspace *w);
0375 int gsl_multifit_robust_maxiter(const size_t maxiter,
0376                                 gsl_multifit_robust_workspace *w);
0377 const char *gsl_multifit_robust_name(const gsl_multifit_robust_workspace *w);
0378 gsl_multifit_robust_stats gsl_multifit_robust_statistics(const gsl_multifit_robust_workspace *w);
0379 int gsl_multifit_robust_weights(const gsl_vector *r, gsl_vector *wts,
0380                                 gsl_multifit_robust_workspace *w);
0381 int gsl_multifit_robust(const gsl_matrix * X, const gsl_vector * y,
0382                         gsl_vector * c, gsl_matrix *cov,
0383                         gsl_multifit_robust_workspace *w);
0384 int gsl_multifit_robust_est(const gsl_vector * x, const gsl_vector * c,
0385                             const gsl_matrix * cov, double *y, double *y_err);
0386 int gsl_multifit_robust_residuals(const gsl_matrix * X,
0387                                   const gsl_vector * y,
0388                                   const gsl_vector * c, gsl_vector * r,
0389                                   gsl_multifit_robust_workspace * w);
0390
0391 __END_DECLS
0392
0393 #endif /* __GSL_MULTIFIT_H__ */