## File: Fitting-robust-linear-regression-example.html

package info (click to toggle)
gsl-ref-html 2.3-1
• area: non-free
• in suites: bullseye, buster, sid
• size: 6,876 kB
• ctags: 4,574
• sloc: makefile: 35
 file content (217 lines) | stat: -rw-r--r-- 7,968 bytes parent folder | download
 `123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217` `````` GNU Scientific Library – Reference Manual: Fitting robust linear regression example

38.8.5 Robust Linear Regression Example

The next program demonstrates the advantage of robust least squares on a dataset with outliers. The program generates linear (x,y) data pairs on the line y = 1.45 x + 3.88, adds some random noise, and inserts 3 outliers into the dataset. Both the robust and ordinary least squares (OLS) coefficients are computed for comparison.

#include <stdio.h> #include <gsl/gsl_multifit.h> #include <gsl/gsl_randist.h>  int dofit(const gsl_multifit_robust_type *T,       const gsl_matrix *X, const gsl_vector *y,       gsl_vector *c, gsl_matrix *cov) {   int s;   gsl_multifit_robust_workspace * work      = gsl_multifit_robust_alloc (T, X->size1, X->size2);    s = gsl_multifit_robust (X, y, c, cov, work);   gsl_multifit_robust_free (work);    return s; }  int main (int argc, char **argv) {   size_t i;   size_t n;   const size_t p = 2; /* linear fit */   gsl_matrix *X, *cov;   gsl_vector *x, *y, *c, *c_ols;   const double a = 1.45; /* slope */   const double b = 3.88; /* intercept */   gsl_rng *r;    if (argc != 2)     {       fprintf (stderr,"usage: robfit n\n");       exit (-1);     }    n = atoi (argv[1]);    X = gsl_matrix_alloc (n, p);   x = gsl_vector_alloc (n);   y = gsl_vector_alloc (n);    c = gsl_vector_alloc (p);   c_ols = gsl_vector_alloc (p);   cov = gsl_matrix_alloc (p, p);    r = gsl_rng_alloc(gsl_rng_default);    /* generate linear dataset */   for (i = 0; i < n - 3; i++)     {       double dx = 10.0 / (n - 1.0);       double ei = gsl_rng_uniform(r);       double xi = -5.0 + i * dx;       double yi = a * xi + b;        gsl_vector_set (x, i, xi);       gsl_vector_set (y, i, yi + ei);     }    /* add a few outliers */   gsl_vector_set(x, n - 3, 4.7);   gsl_vector_set(y, n - 3, -8.3);    gsl_vector_set(x, n - 2, 3.5);   gsl_vector_set(y, n - 2, -6.7);    gsl_vector_set(x, n - 1, 4.1);   gsl_vector_set(y, n - 1, -6.0);    /* construct design matrix X for linear fit */   for (i = 0; i < n; ++i)     {       double xi = gsl_vector_get(x, i);        gsl_matrix_set (X, i, 0, 1.0);       gsl_matrix_set (X, i, 1, xi);     }    /* perform robust and OLS fit */   dofit(gsl_multifit_robust_ols, X, y, c_ols, cov);   dofit(gsl_multifit_robust_bisquare, X, y, c, cov);    /* output data and model */   for (i = 0; i < n; ++i)     {       double xi = gsl_vector_get(x, i);       double yi = gsl_vector_get(y, i);       gsl_vector_view v = gsl_matrix_row(X, i);       double y_ols, y_rob, y_err;        gsl_multifit_robust_est(&v.vector, c, cov, &y_rob, &y_err);       gsl_multifit_robust_est(&v.vector, c_ols, cov, &y_ols, &y_err);        printf("%g %g %g %g\n", xi, yi, y_rob, y_ols);     }  #define C(i) (gsl_vector_get(c,(i))) #define COV(i,j) (gsl_matrix_get(cov,(i),(j)))    {     printf ("# best fit: Y = %g + %g X\n",              C(0), C(1));      printf ("# covariance matrix:\n");     printf ("# [ %+.5e, %+.5e\n",                COV(0,0), COV(0,1));     printf ("#   %+.5e, %+.5e\n",                 COV(1,0), COV(1,1));   }    gsl_matrix_free (X);   gsl_vector_free (x);   gsl_vector_free (y);   gsl_vector_free (c);   gsl_vector_free (c_ols);   gsl_matrix_free (cov);   gsl_rng_free(r);    return 0; }

The output from the program is shown in the following plot.

``````