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

   1 size_t pontius_n = 40;
   2 size_t pontius_p = 3;
   3
   4 double pontius_x[] = { 150000, 300000, 450000, 600000, 750000, 900000,
   5 1050000, 1200000, 1350000, 1500000, 1650000, 1800000, 1950000, 2100000,
   6 2250000, 2400000, 2550000, 2700000, 2850000, 3000000, 150000, 300000,
   7 450000, 600000, 750000, 900000, 1050000, 1200000, 1350000, 1500000,
   8 1650000, 1800000, 1950000, 2100000, 2250000, 2400000, 2550000, 2700000,
   9 2850000, 3000000 };
  10
  11 double pontius_y[] = { .11019, .21956, .32949, .43899, .54803, .65694,
  12 .76562, .87487, .98292, 1.09146, 1.20001, 1.30822, 1.41599, 1.52399,
  13 1.63194, 1.73947, 1.84646, 1.95392, 2.06128, 2.16844, .11052, .22018,
  14 .32939, .43886, .54798, .65739, .76596, .87474, .98300, 1.09150,
  15 1.20004, 1.30818, 1.41613, 1.52408, 1.63159, 1.73965, 1.84696,
  16 1.95445, 2.06177, 2.16829 };
  17
  18 void
  19 test_pontius ()
  20 {
  21   size_t i, j;
  22   {
  23     gsl_multifit_linear_workspace * work =
  24       gsl_multifit_linear_alloc (pontius_n, pontius_p);
  25
  26     gsl_matrix * X = gsl_matrix_alloc (pontius_n, pontius_p);
  27     gsl_vector_view y = gsl_vector_view_array (pontius_y, pontius_n);
  28     gsl_vector * c = gsl_vector_alloc (pontius_p);
  29     gsl_vector * r = gsl_vector_alloc (pontius_n);
  30     gsl_matrix * cov = gsl_matrix_alloc (pontius_p, pontius_p);
  31     gsl_vector_view diag;
  32
  33     double chisq;
  34
  35     double expected_c[3] = { 0.673565789473684E-03,
  36                              0.732059160401003E-06,
  37                             -0.316081871345029E-14};
  38
  39     double expected_sd[3] = { 0.107938612033077E-03,
  40                               0.157817399981659E-09,
  41                               0.486652849992036E-16 };
  42
  43     double expected_chisq = 0.155761768796992E-05;
  44
  45     for (i = 0 ; i < pontius_n; i++)
  46       {
  47         for (j = 0; j < pontius_p; j++)
  48           {
  49             gsl_matrix_set(X, i, j, pow(pontius_x[i], j));
  50           }
  51       }
  52
  53     gsl_multifit_linear (X, &y.vector, c, cov, &chisq, work);
  54
  55     gsl_test_rel (gsl_vector_get(c,0), expected_c[0], 1e-10, "pontius gsl_fit_multilinear c0") ;
  56     gsl_test_rel (gsl_vector_get(c,1), expected_c[1], 1e-10, "pontius gsl_fit_multilinear c1") ;
  57     gsl_test_rel (gsl_vector_get(c,2), expected_c[2], 1e-10, "pontius gsl_fit_multilinear c2") ;
  58
  59     diag = gsl_matrix_diagonal (cov);
  60
  61     gsl_test_rel (gsl_vector_get(&diag.vector,0), pow(expected_sd[0],2.0), 1e-10, "pontius gsl_fit_multilinear cov00") ;
  62     gsl_test_rel (gsl_vector_get(&diag.vector,1), pow(expected_sd[1],2.0), 1e-10, "pontius gsl_fit_multilinear cov11") ;
  63     gsl_test_rel (gsl_vector_get(&diag.vector,2), pow(expected_sd[2],2.0), 1e-10, "pontius gsl_fit_multilinear cov22") ;
  64
  65     gsl_test_rel (chisq, expected_chisq, 1e-10, "pontius gsl_fit_multilinear chisq") ;
  66
  67     gsl_multifit_linear_residuals(X, &y.vector, c, r);
  68     gsl_blas_ddot(r, r, &chisq);
  69     gsl_test_rel (chisq, expected_chisq, 1e-10, "pontius gsl_fit_multilinear residuals") ;
  70
  71     gsl_vector_free(c);
  72     gsl_vector_free(r);
  73     gsl_matrix_free(cov);
  74     gsl_matrix_free(X);
  75     gsl_multifit_linear_free (work);
  76   }
  77
  78
  79   {
  80     gsl_multifit_linear_workspace * work =
  81       gsl_multifit_linear_alloc (pontius_n, pontius_p);
  82
  83     gsl_matrix * X = gsl_matrix_alloc (pontius_n, pontius_p);
  84     gsl_vector_view y = gsl_vector_view_array (pontius_y, pontius_n);
  85     gsl_vector * w = gsl_vector_alloc (pontius_n);
  86     gsl_vector * c = gsl_vector_alloc (pontius_p);
  87     gsl_vector * r = gsl_vector_alloc (pontius_n);
  88     gsl_matrix * cov = gsl_matrix_alloc (pontius_p, pontius_p);
  89
  90     double chisq;
  91
  92     double expected_c[3] = {  0.673565789473684E-03,
  93                                0.732059160401003E-06,
  94                                -0.316081871345029E-14};
  95
  96     double expected_chisq = 0.155761768796992E-05;
  97
  98     double expected_cov[3][3] ={
  99       {2.76754385964916e-01 , -3.59649122807024e-07,   9.74658869395731e-14},
 100       {-3.59649122807024e-07,   5.91630591630603e-13,  -1.77210703526497e-19},
 101       {9.74658869395731e-14,  -1.77210703526497e-19,   5.62573661988878e-26} };
 102
 103
 104     for (i = 0 ; i < pontius_n; i++)
 105       {
 106         for (j = 0; j < pontius_p; j++)
 107           {
 108             gsl_matrix_set(X, i, j, pow(pontius_x[i], j));
 109           }
 110       }
 111
 112     gsl_vector_set_all (w, 1.0);
 113
 114     gsl_multifit_wlinear (X, w, &y.vector, c, cov, &chisq, work);
 115
 116     gsl_test_rel (gsl_vector_get(c,0), expected_c[0], 1e-10, "pontius gsl_fit_multilinear c0") ;
 117     gsl_test_rel (gsl_vector_get(c,1), expected_c[1], 1e-10, "pontius gsl_fit_multilinear c1") ;
 118     gsl_test_rel (gsl_vector_get(c,2), expected_c[2], 1e-10, "pontius gsl_fit_multilinear c2") ;
 119
 120
 121     for (i = 0; i < pontius_p; i++)
 122       {
 123         for (j = 0; j < pontius_p; j++)
 124           {
 125             gsl_test_rel (gsl_matrix_get(cov,i,j), expected_cov[i][j], 1e-10,
 126                           "pontius gsl_fit_wmultilinear cov(%d,%d)", i, j) ;
 127           }
 128       }
 129
 130     gsl_test_rel (chisq, expected_chisq, 1e-10, "pontius gsl_fit_wmultilinear chisq") ;
 131
 132     gsl_multifit_linear_residuals(X, &y.vector, c, r);
 133     gsl_blas_ddot(r, r, &chisq);
 134     gsl_test_rel (chisq, expected_chisq, 1e-10, "pontius gsl_fit_wmultilinear residuals") ;
 135
 136     gsl_vector_free(w);
 137     gsl_vector_free(c);
 138     gsl_vector_free(r);
 139     gsl_matrix_free(cov);
 140     gsl_matrix_free(X);
 141     gsl_multifit_linear_free (work);
 142   }
 143 }