source: trunk/test/regression_test.cc @ 730

Last change on this file since 730 was 730, checked in by Peter, 16 years ago

fixes #167 and #160

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 10.7 KB
Line 
1// $Id: regression_test.cc 730 2007-01-06 11:02:21Z peter $
2
3/*
4  Copyright (C) The authors contributing to this file.
5
6  This file is part of the yat library, http://lev.thep.lu.se/trac/yat
7
8  The yat library is free software; you can redistribute it and/or
9  modify it under the terms of the GNU General Public License as
10  published by the Free Software Foundation; either version 2 of the
11  License, or (at your option) any later version.
12
13  The yat library is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16  General Public License for more details.
17
18  You should have received a copy of the GNU General Public License
19  along with this program; if not, write to the Free Software
20  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
21  02111-1307, USA.
22*/
23
24#include "yat/regression/KernelBox.h"
25#include "yat/regression/Linear.h"
26#include "yat/regression/LinearWeighted.h"
27#include "yat/regression/Local.h"
28#include "yat/regression/Naive.h"
29#include "yat/regression/NaiveWeighted.h"
30#include "yat/regression/Polynomial.h"
31#include "yat/regression/PolynomialWeighted.h"
32#include "yat/utility/matrix.h"
33#include "yat/utility/vector.h"
34
35#include <cmath>
36
37#include <fstream>
38#include <iostream>
39
40
41using namespace theplu::yat;
42
43bool equal(regression::OneDimensional&, regression::OneDimensionalWeighted&, 
44           std::ostream*); 
45
46bool unity_weights(regression::OneDimensional&, 
47                   regression::OneDimensionalWeighted&, 
48                   const utility::vector&, const utility::vector&,
49                   std::ostream*); 
50
51bool rescale_weights(regression::OneDimensionalWeighted&, 
52                     const utility::vector&, const utility::vector&,
53                     std::ostream*); 
54
55bool zero_weights(regression::OneDimensionalWeighted&, 
56                  const utility::vector&, const utility::vector&,
57                  std::ostream*); 
58
59
60bool Local_test(regression::OneDimensionalWeighted&, 
61                regression::Kernel&);
62
63int main(const int argc,const char* argv[])
64{
65  std::ostream* error;
66  if (argc>1 && argv[1]==std::string("-v"))
67    error = &std::cerr;
68  else {
69    error = new std::ofstream("/dev/null");
70    if (argc>1)
71      std::cout << "regression_test -v : for printing extra information\n";
72  }
73  *error << "  testing regression" << std::endl;
74  bool ok = true;
75
76  // test data for Linear and Naive (Weighted and non-weighted)
77  utility::vector x(4); x(0)=1970; x(1)=1980; x(2)=1990; x(3)=2000;
78  utility::vector y(4); y(0)=12;   y(1)=11;   y(2)=14;   y(3)=13;
79  utility::vector w(4); w(0)=0.1;  w(1)=0.2;  w(2)=0.3;  w(3)=0.4;
80
81  // Comparing linear and polynomial(1)
82  regression::Linear linear;
83  linear.fit(x,y);
84  regression::Polynomial polynomial(1);
85  polynomial.fit(x,y);
86  if ( fabs(linear.beta()-polynomial.fit_parameters()(1))>0.0001 ){
87    *error << "error: beta and fit_parameters(1) not equal" << std::endl;
88    *error << "       beta = " << linear.beta() << std::endl;
89    *error << "       fit_parameters(1) = " 
90           << polynomial.fit_parameters()(1) << std::endl;
91    ok = false;
92  }
93  if ( fabs(polynomial.fit_parameters()(0)-linear.alpha()+
94            linear.beta()*1985)>0.0001){
95    *error << "error: fit_parameters(0) = " 
96           << polynomial.fit_parameters()(0)<< std::endl;
97    *error << "error: alpha-beta*m_x = " 
98           << linear.alpha()-linear.beta()*1985 << std::endl;
99    ok = false;
100  }
101  if ( fabs(polynomial.chisq()-linear.chisq())>0.0001){
102    *error << "error: chisq not same in linear and polynomial(1)" 
103           << std::endl;
104    ok = false;
105  }
106  if ( fabs(polynomial.predict(1.0)-linear.predict(1.0))>0.0001){
107    *error << "error: predict not same in linear and polynomial(1)" 
108           << std::endl;
109    ok = false;
110  }
111  if ( fabs(polynomial.standard_error2(1985)-linear.standard_error2(1985))
112       >0.0001){
113    *error << "error: standard_error not same in linear and polynomial(1)" 
114           << "\n  polynomial: " << polynomial.standard_error2(1985)
115           << "\n  linear: " << linear.standard_error2(1985)
116           << std::endl;
117    ok = false;
118  }
119
120  *error << "  testing regression::LinearWeighted" << std::endl;
121  regression::LinearWeighted linear_w;
122  ok = equal(linear, linear_w, error) && ok;
123  linear_w.fit(x,y,w);
124  double y_predicted = linear_w.predict(1990);
125  double y_predicted_err = linear_w.prediction_error2(1990);
126  if (fabs(y_predicted-12.8)>0.001){
127    *error << "error: cannot reproduce fit." << std::endl;
128    *error << "predicted value: " << y_predicted << " expected 12.8" 
129           << std::endl;
130    ok=false;
131  }
132
133  // testing regression::NaiveWeighted
134  *error << "  testing regression::NaiveWeighted" << std::endl;
135  regression::NaiveWeighted naive_w;
136  regression::Naive naive;
137  ok = equal(naive, naive_w, error) && ok;
138  naive_w.fit(x,y,w);
139
140  y_predicted=naive_w.predict(0.0);
141  y_predicted_err=naive_w.prediction_error2(0.0);
142  if (y_predicted!=(0.1*12+0.2*11+0.3*14+0.4*13)) {
143    *error << "regression_NaiveWeighted: cannot reproduce fit." << std::endl;
144    *error << "returned value: " << y_predicted << std::endl;
145    *error << "expected: " << 0.1*12+0.2*11+0.3*14+0.4*13 << std::endl;
146    ok=false;
147  }
148
149  // testing regression::Local
150  *error << "  testing regression::Local" << std::endl;
151  regression::KernelBox kb;
152  regression::LinearWeighted rl;
153  if (!Local_test(rl,kb)) {
154    *error << "regression_Local: Linear cannot reproduce fit." << std::endl;
155    ok=false;
156  }
157  regression::NaiveWeighted rn;
158  if (!Local_test(rn,kb)) {
159    *error << "regression_Local: Naive cannot reproduce fit." << std::endl;
160    ok=false;
161  }
162
163  // testing regression::Polynomial
164  *error << "  testing regression::Polynomial" << std::endl;
165  {
166    std::ifstream s("data/regression_gauss.data");
167    utility::matrix data(s);
168    utility::vector x(data.rows());
169    utility::vector ln_y(data.rows());
170    for (size_t i=0; i<data.rows(); ++i) {
171      x(i)=data(i,0);
172      ln_y(i)=log(data(i,1));
173    }
174
175    regression::Polynomial polynomialfit(2);
176    polynomialfit.fit(x,ln_y);
177    utility::vector fit=polynomialfit.fit_parameters();
178    if (fabs(fit[0]-1.012229646706 + fit[1]-0.012561322528 +
179             fit[2]+1.159674470130)>1e-11) {  // Jari, fix number!
180      *error << "regression_Polynomial: cannot reproduce fit." << std::endl;
181      ok=false;
182    }
183  }
184
185  *error << "  testing regression::PolynomialWeighted" << std::endl;
186  regression::PolynomialWeighted pol_weighted(2);
187
188  if (error!=&std::cerr)
189    delete error;
190
191  return (ok ? 0 : -1);
192}
193
194
195bool equal(regression::OneDimensional& r, 
196           regression::OneDimensionalWeighted& wr, 
197           std::ostream* error)
198{
199  bool ok=true;
200  utility::vector x(5); x(0)=1970; x(1)=1980; x(2)=1990; x(3)=2000; x(4)=2010;
201  utility::vector y(5); y(0)=12;   y(1)=11;   y(2)=14;   y(3)=13;   y(4)=15;
202
203  ok = unity_weights(r, wr, x, y, error) && ok;
204  ok = rescale_weights(wr, x, y, error) && ok;
205  ok = zero_weights(wr, x, y, error) && ok;
206  return ok;
207}
208
209
210bool unity_weights(regression::OneDimensional& r, 
211                   regression::OneDimensionalWeighted& rw, 
212                   const utility::vector& x, const utility::vector& y,
213                   std::ostream* error)
214{
215  *error << "    testing unity weights equal to non-weighted version.\n"; 
216  bool ok=true;
217  utility::vector w(x.size(), 1.0);
218  r.fit(x,y);
219  rw.fit(x,y,w);
220  if (r.predict(2000) != rw.predict(2000)){
221    ok = false;
222    *error << "Error: predict not equal" << std::endl;
223  }
224  if (fabs(r.prediction_error2(2000)-rw.prediction_error2(2000))>10e-7){
225    ok = false;
226    *error << "Error: prediction_error2 not equal non-weighted version.\n" 
227           << "       weighted: " << rw.prediction_error2(2000) << "\n"
228           << "   non-weighted: " << r.prediction_error2(2000)
229           << std::endl;
230  }
231  if (fabs(r.s2()-rw.s2())>10E-7){
232    ok = false;
233    *error << "Error: r2 not equal non-weighted version." << std::endl;
234    *error << "weighted r2 = " << rw.r2() << std::endl;
235    *error << "non-weighted r2 = " << r.r2() << std::endl;
236  }
237  if (fabs(r.s2()-rw.s2())>10E-7){
238    ok = false;
239    *error << "Error: s2 not equal non-weighted version." << std::endl;
240    *error << "weighted s2 = " << rw.s2() << std::endl;
241    *error << "non-weighted s2 = " << r.s2() << std::endl;
242  }
243  if (fabs(r.standard_error2(2000)-rw.standard_error2(2000))>10e-7){
244    ok = false;
245    *error << "Error: standard_error not equal non-weighted version." 
246           << std::endl;
247  }
248  return ok;
249} 
250
251
252bool rescale_weights(regression::OneDimensionalWeighted& wr, 
253                     const utility::vector& x, const utility::vector& y,
254                     std::ostream* error)
255{
256  *error << "    testing rescaling weights.\n"; 
257  bool ok = true;
258  utility::vector w(5);  w(0)=1.0;  w(1)=1.0;  w(2)=0.5;  w(3)=0.2;  w(4)=0.2;
259  wr.fit(x,y,w);
260  double predict = wr.predict(2000);
261  double prediction_error2 = wr.prediction_error2(2000);
262  double r2 = wr.r2();
263  double s2 = wr.s2();
264  double standard_error2 = wr.standard_error2(2000);
265
266  w.scale(2);
267  wr.fit(x,y,w);
268  if (wr.predict(2000) != predict){
269    ok = false;
270    *error << "Error: predict not equal after rescaling.\n";
271  }
272  if (wr.prediction_error2(2000,2) != prediction_error2){
273    ok = false;
274    *error << "Error: prediction_error2 not equal after rescaling.\n";
275  }
276  if (wr.r2() != r2){
277    ok = false;
278    *error << "Error: r2 not equal after rescaling.\n";
279  }
280  if (wr.s2(2) != s2){
281    ok = false;
282    *error << "Error: s2 not equal after rescaling.\n";
283    *error << "       s2 = " << s2 << " and after doubling weights.\n";
284    *error << "       s2 = " << wr.s2(2) << "\n";
285  }
286  if (wr.standard_error2(2000) != standard_error2){
287    ok = false;
288    *error << "Error: standard_error2 not equal after rescaling.\n";
289  }
290  return ok; 
291}
292
293
294bool zero_weights(regression::OneDimensionalWeighted& wr, 
295                  const utility::vector& x, const utility::vector& y,
296                  std::ostream* error)
297{
298  *error << "    testing zero weights equal to missing value.\n"; 
299  bool ok = true;
300  utility::vector w(5);  w(0)=1.0;  w(1)=1.0;  w(2)=0.5;  w(3)=0.2;  w(4)=0;
301  wr.fit(x,y,w);
302  double predict = wr.predict(2000);
303  double prediction_error2 = wr.prediction_error2(2000);
304  double r2 = wr.r2();
305  double s2 = wr.s2();
306  double standard_error2 = wr.standard_error2(2000);
307
308  utility::vector x2(4);
309  utility::vector y2(4);
310  utility::vector w2(4);
311  for (size_t i=0; i<4; ++i){
312    x2(i) = x(i);
313    y2(i) = y(i);
314    w2(i) = w(i);
315  }
316
317  wr.fit(x2,y2,w2);
318  if (wr.predict(2000) != predict){
319    ok = false;
320    *error << "Error: predict not equal.\n";
321  }
322  if (wr.prediction_error2(2000) != prediction_error2){
323    ok = false;
324    *error << "Error: prediction_error2 not equal.\n";
325  }
326  if (wr.r2() != r2){
327    ok = false;
328    *error << "Error: r2 not equal.\n";
329  }
330  if (wr.s2() != s2){
331    ok = false;
332    *error << "Error: s2 not equal.\n";
333  }
334  if (wr.standard_error2(2000) != standard_error2){
335    ok = false;
336    *error << "Error: standard_error2 not equal.\n";
337  }
338  return ok; 
339}
340
341
342bool Local_test(regression::OneDimensionalWeighted& r, 
343                regression::Kernel& k)
344{
345  regression::Local rl(r,k);
346  for (size_t i=0; i<1000; i++){
347    rl.add(i, 10);
348  }
349
350  rl.fit(10, 100);
351
352  utility::vector y = rl.y_predicted();
353  for (size_t i=0; i<y.size(); i++) 
354    if (y(i)!=10.0){
355      return false;
356    }
357  return true;
358}
Note: See TracBrowser for help on using the repository browser.