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

   1 size_t longley_n = 16;
   2 size_t longley_p = 7;
   3
   4 double longley_x [] = {
   5   1,  83.0,   234289,   2356,     1590,    107608,  1947,
   6   1,  88.5,   259426,   2325,     1456,    108632,  1948,
   7   1,  88.2,   258054,   3682,     1616,    109773,  1949,
   8   1,  89.5,   284599,   3351,     1650,    110929,  1950,
   9   1,  96.2,   328975,   2099,     3099,    112075,  1951,
  10   1,  98.1,   346999,   1932,     3594,    113270,  1952,
  11   1,  99.0,   365385,   1870,     3547,    115094,  1953,
  12   1, 100.0,   363112,   3578,     3350,    116219,  1954,
  13   1, 101.2,   397469,   2904,     3048,    117388,  1955,
  14   1, 104.6,   419180,   2822,     2857,    118734,  1956,
  15   1, 108.4,   442769,   2936,     2798,    120445,  1957,
  16   1, 110.8,   444546,   4681,     2637,    121950,  1958,
  17   1, 112.6,   482704,   3813,     2552,    123366,  1959,
  18   1, 114.2,   502601,   3931,     2514,    125368,  1960,
  19   1, 115.7,   518173,   4806,     2572,    127852,  1961,
  20   1, 116.9,   554894,   4007,     2827,    130081,  1962 } ;
  21
  22 double longley_y[] = {60323, 61122, 60171, 61187, 63221, 63639, 64989, 63761,
  23                        66019, 67857, 68169, 66513, 68655, 69564, 69331, 70551};
  24
  25
  26 void
  27 test_longley ()
  28 {
  29   size_t i, j;
  30   {
  31     gsl_multifit_linear_workspace * work =
  32       gsl_multifit_linear_alloc (longley_n, longley_p);
  33
  34     gsl_matrix_view X = gsl_matrix_view_array (longley_x, longley_n, longley_p);
  35     gsl_vector_view y = gsl_vector_view_array (longley_y, longley_n);
  36     gsl_vector * c = gsl_vector_alloc (longley_p);
  37     gsl_vector * r = gsl_vector_alloc (longley_n);
  38     gsl_matrix * cov = gsl_matrix_alloc (longley_p, longley_p);
  39     gsl_vector_view diag;
  40
  41     double chisq;
  42
  43     double expected_c[7] = {  -3482258.63459582,
  44                               15.0618722713733,
  45                               -0.358191792925910E-01,
  46                               -2.02022980381683,
  47                               -1.03322686717359,
  48                               -0.511041056535807E-01,
  49                               1829.15146461355 };
  50
  51     double expected_sd[7]  = {  890420.383607373,
  52                                 84.9149257747669,
  53                                 0.334910077722432E-01,
  54                                 0.488399681651699,
  55                                 0.214274163161675,
  56                                 0.226073200069370,
  57                                 455.478499142212 } ;
  58
  59     double expected_chisq = 836424.055505915;
  60
  61     gsl_multifit_linear (&X.matrix, &y.vector, c, cov, &chisq, work);
  62
  63     gsl_test_rel (gsl_vector_get(c,0), expected_c[0], 1e-10, "longley gsl_fit_multilinear c0") ;
  64     gsl_test_rel (gsl_vector_get(c,1), expected_c[1], 1e-10, "longley gsl_fit_multilinear c1") ;
  65     gsl_test_rel (gsl_vector_get(c,2), expected_c[2], 1e-10, "longley gsl_fit_multilinear c2") ;
  66     gsl_test_rel (gsl_vector_get(c,3), expected_c[3], 1e-10, "longley gsl_fit_multilinear c3") ;
  67     gsl_test_rel (gsl_vector_get(c,4), expected_c[4], 1e-10, "longley gsl_fit_multilinear c4") ;
  68     gsl_test_rel (gsl_vector_get(c,5), expected_c[5], 1e-10, "longley gsl_fit_multilinear c5") ;
  69     gsl_test_rel (gsl_vector_get(c,6), expected_c[6], 1e-10, "longley gsl_fit_multilinear c6") ;
  70
  71     diag = gsl_matrix_diagonal (cov);
  72
  73     gsl_test_rel (gsl_vector_get(&diag.vector,0), pow(expected_sd[0],2.0), 1e-10, "longley gsl_fit_multilinear cov00") ;
  74     gsl_test_rel (gsl_vector_get(&diag.vector,1), pow(expected_sd[1],2.0), 1e-10, "longley gsl_fit_multilinear cov11") ;
  75     gsl_test_rel (gsl_vector_get(&diag.vector,2), pow(expected_sd[2],2.0), 1e-10, "longley gsl_fit_multilinear cov22") ;
  76     gsl_test_rel (gsl_vector_get(&diag.vector,3), pow(expected_sd[3],2.0), 1e-10, "longley gsl_fit_multilinear cov33") ;
  77     gsl_test_rel (gsl_vector_get(&diag.vector,4), pow(expected_sd[4],2.0), 1e-10, "longley gsl_fit_multilinear cov44") ;
  78     gsl_test_rel (gsl_vector_get(&diag.vector,5), pow(expected_sd[5],2.0), 1e-10, "longley gsl_fit_multilinear cov55") ;
  79     gsl_test_rel (gsl_vector_get(&diag.vector,6), pow(expected_sd[6],2.0), 1e-10, "longley gsl_fit_multilinear cov66") ;
  80
  81     gsl_test_rel (chisq, expected_chisq, 1e-10, "longley gsl_fit_multilinear chisq") ;
  82
  83     gsl_multifit_linear_residuals(&X.matrix, &y.vector, c, r);
  84     gsl_blas_ddot(r, r, &chisq);
  85     gsl_test_rel (chisq, expected_chisq, 1e-10, "longley gsl_fit_multilinear residuals") ;
  86
  87     gsl_vector_free(c);
  88     gsl_vector_free(r);
  89     gsl_matrix_free(cov);
  90     gsl_multifit_linear_free (work);
  91   }
  92
  93
  94   {
  95     gsl_multifit_linear_workspace * work =
  96       gsl_multifit_linear_alloc (longley_n, longley_p);
  97
  98     gsl_matrix_view X = gsl_matrix_view_array (longley_x, longley_n, longley_p);
  99     gsl_vector_view y = gsl_vector_view_array (longley_y, longley_n);
 100     gsl_vector * w = gsl_vector_alloc (longley_n);
 101     gsl_vector * c = gsl_vector_alloc (longley_p);
 102     gsl_vector * r = gsl_vector_alloc (longley_n);
 103     gsl_matrix * cov = gsl_matrix_alloc (longley_p, longley_p);
 104
 105     double chisq;
 106
 107     double expected_c[7] = {  -3482258.63459582,
 108                               15.0618722713733,
 109                               -0.358191792925910E-01,
 110                               -2.02022980381683,
 111                               -1.03322686717359,
 112                               -0.511041056535807E-01,
 113                               1829.15146461355 };
 114
 115     double expected_cov[7][7] = { { 8531122.56783558,
 116 -166.727799925578, 0.261873708176346, 3.91188317230983,
 117 1.1285582054705, -0.889550869422687, -4362.58709870581},
 118
 119 {-166.727799925578, 0.0775861253030891, -1.98725210399982e-05,
 120 -0.000247667096727256, -6.82911920718824e-05, 0.000136160797527761,
 121 0.0775255245956248},
 122
 123 {0.261873708176346, -1.98725210399982e-05, 1.20690316701888e-08,
 124 1.66429546772984e-07, 3.61843600487847e-08, -6.78805814483582e-08,
 125 -0.00013158719037715},
 126
 127 {3.91188317230983, -0.000247667096727256, 1.66429546772984e-07,
 128 2.56665052544717e-06, 6.96541409215597e-07, -9.00858307771567e-07,
 129 -0.00197260370663974},
 130
 131 {1.1285582054705, -6.82911920718824e-05, 3.61843600487847e-08,
 132 6.96541409215597e-07, 4.94032602583969e-07, -9.8469143760973e-08,
 133 -0.000576921112208274},
 134
 135 {-0.889550869422687, 0.000136160797527761, -6.78805814483582e-08,
 136 -9.00858307771567e-07, -9.8469143760973e-08, 5.49938542664952e-07,
 137 0.000430074434198215},
 138
 139 {-4362.58709870581, 0.0775255245956248, -0.00013158719037715,
 140 -0.00197260370663974, -0.000576921112208274, 0.000430074434198215,
 141 2.23229587481535 }} ;
 142
 143     double expected_chisq = 836424.055505915;
 144
 145     gsl_vector_set_all (w, 1.0);
 146
 147     gsl_multifit_wlinear (&X.matrix, w, &y.vector, c, cov, &chisq, work);
 148
 149     gsl_test_rel (gsl_vector_get(c,0), expected_c[0], 1e-10, "longley gsl_fit_wmultilinear c0") ;
 150     gsl_test_rel (gsl_vector_get(c,1), expected_c[1], 1e-10, "longley gsl_fit_wmultilinear c1") ;
 151     gsl_test_rel (gsl_vector_get(c,2), expected_c[2], 1e-10, "longley gsl_fit_wmultilinear c2") ;
 152     gsl_test_rel (gsl_vector_get(c,3), expected_c[3], 1e-10, "longley gsl_fit_wmultilinear c3") ;
 153     gsl_test_rel (gsl_vector_get(c,4), expected_c[4], 1e-10, "longley gsl_fit_wmultilinear c4") ;
 154     gsl_test_rel (gsl_vector_get(c,5), expected_c[5], 1e-10, "longley gsl_fit_wmultilinear c5") ;
 155     gsl_test_rel (gsl_vector_get(c,6), expected_c[6], 1e-10, "longley gsl_fit_wmultilinear c6") ;
 156
 157     for (i = 0; i < longley_p; i++)
 158       {
 159         for (j = 0; j < longley_p; j++)
 160           {
 161             gsl_test_rel (gsl_matrix_get(cov,i,j), expected_cov[i][j], 1e-7,
 162                           "longley gsl_fit_wmultilinear cov(%d,%d)", i, j) ;
 163           }
 164       }
 165
 166     gsl_test_rel (chisq, expected_chisq, 1e-10, "longley gsl_fit_wmultilinear chisq") ;
 167
 168     gsl_multifit_linear_residuals(&X.matrix, &y.vector, c, r);
 169     gsl_blas_ddot(r, r, &chisq);
 170     gsl_test_rel (chisq, expected_chisq, 1e-10, "longley gsl_fit_wmultilinear residuals") ;
 171
 172     gsl_vector_free(w);
 173     gsl_vector_free(c);
 174     gsl_vector_free(r);
 175     gsl_matrix_free(cov);
 176     gsl_multifit_linear_free (work);
 177   }
 178 }