source: trunk/test/regression_test.cc @ 731

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

added test for multidimensional weighted and straight version

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 12.0 KB
Line 
1// $Id: regression_test.cc 731 2007-01-06 16:06:19Z 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 multidim(std::ostream* error);
47
48bool unity_weights(regression::OneDimensional&, 
49                   regression::OneDimensionalWeighted&, 
50                   const utility::vector&, const utility::vector&,
51                   std::ostream*); 
52
53bool rescale_weights(regression::OneDimensionalWeighted&, 
54                     const utility::vector&, const utility::vector&,
55                     std::ostream*); 
56
57bool zero_weights(regression::OneDimensionalWeighted&, 
58                  const utility::vector&, const utility::vector&,
59                  std::ostream*); 
60
61
62bool Local_test(regression::OneDimensionalWeighted&, 
63                regression::Kernel&);
64
65int main(const int argc,const char* argv[])
66{
67  std::ostream* error;
68  if (argc>1 && argv[1]==std::string("-v"))
69    error = &std::cerr;
70  else {
71    error = new std::ofstream("/dev/null");
72    if (argc>1)
73      std::cout << "regression_test -v : for printing extra information\n";
74  }
75  *error << "  testing regression" << std::endl;
76  bool ok = true;
77
78  // test data for Linear and Naive (Weighted and non-weighted)
79  utility::vector x(4); x(0)=1970; x(1)=1980; x(2)=1990; x(3)=2000;
80  utility::vector y(4); y(0)=12;   y(1)=11;   y(2)=14;   y(3)=13;
81  utility::vector w(4); w(0)=0.1;  w(1)=0.2;  w(2)=0.3;  w(3)=0.4;
82
83  // Comparing linear and polynomial(1)
84  regression::Linear linear;
85  linear.fit(x,y);
86  regression::Polynomial polynomial(1);
87  polynomial.fit(x,y);
88  if ( fabs(linear.beta()-polynomial.fit_parameters()(1))>0.0001 ){
89    *error << "error: beta and fit_parameters(1) not equal" << std::endl;
90    *error << "       beta = " << linear.beta() << std::endl;
91    *error << "       fit_parameters(1) = " 
92           << polynomial.fit_parameters()(1) << std::endl;
93    ok = false;
94  }
95  if ( fabs(polynomial.fit_parameters()(0)-linear.alpha()+
96            linear.beta()*1985)>0.0001){
97    *error << "error: fit_parameters(0) = " 
98           << polynomial.fit_parameters()(0)<< std::endl;
99    *error << "error: alpha-beta*m_x = " 
100           << linear.alpha()-linear.beta()*1985 << std::endl;
101    ok = false;
102  }
103  if ( fabs(polynomial.chisq()-linear.chisq())>0.0001){
104    *error << "error: chisq not same in linear and polynomial(1)" 
105           << std::endl;
106    ok = false;
107  }
108  if ( fabs(polynomial.predict(1.0)-linear.predict(1.0))>0.0001){
109    *error << "error: predict not same in linear and polynomial(1)" 
110           << std::endl;
111    ok = false;
112  }
113  if ( fabs(polynomial.standard_error2(1985)-linear.standard_error2(1985))
114       >0.0001){
115    *error << "error: standard_error not same in linear and polynomial(1)" 
116           << "\n  polynomial: " << polynomial.standard_error2(1985)
117           << "\n  linear: " << linear.standard_error2(1985)
118           << std::endl;
119    ok = false;
120  }
121
122  *error << "  testing regression::LinearWeighted" << std::endl;
123  regression::LinearWeighted linear_w;
124  ok = equal(linear, linear_w, error) && ok;
125  linear_w.fit(x,y,w);
126  double y_predicted = linear_w.predict(1990);
127  double y_predicted_err = linear_w.prediction_error2(1990);
128  if (fabs(y_predicted-12.8)>0.001){
129    *error << "error: cannot reproduce fit." << std::endl;
130    *error << "predicted value: " << y_predicted << " expected 12.8" 
131           << std::endl;
132    ok=false;
133  }
134
135  // testing regression::NaiveWeighted
136  *error << "  testing regression::NaiveWeighted" << std::endl;
137  regression::NaiveWeighted naive_w;
138  regression::Naive naive;
139  ok = equal(naive, naive_w, error) && ok;
140  naive_w.fit(x,y,w);
141
142  y_predicted=naive_w.predict(0.0);
143  y_predicted_err=naive_w.prediction_error2(0.0);
144  if (y_predicted!=(0.1*12+0.2*11+0.3*14+0.4*13)) {
145    *error << "regression_NaiveWeighted: cannot reproduce fit." << std::endl;
146    *error << "returned value: " << y_predicted << std::endl;
147    *error << "expected: " << 0.1*12+0.2*11+0.3*14+0.4*13 << std::endl;
148    ok=false;
149  }
150
151  // testing regression::Local
152  *error << "  testing regression::Local" << std::endl;
153  regression::KernelBox kb;
154  regression::LinearWeighted rl;
155  if (!Local_test(rl,kb)) {
156    *error << "regression_Local: Linear cannot reproduce fit." << std::endl;
157    ok=false;
158  }
159  regression::NaiveWeighted rn;
160  if (!Local_test(rn,kb)) {
161    *error << "regression_Local: Naive cannot reproduce fit." << std::endl;
162    ok=false;
163  }
164
165  // testing regression::Polynomial
166  *error << "  testing regression::Polynomial" << std::endl;
167  {
168    std::ifstream s("data/regression_gauss.data");
169    utility::matrix data(s);
170    utility::vector x(data.rows());
171    utility::vector ln_y(data.rows());
172    for (size_t i=0; i<data.rows(); ++i) {
173      x(i)=data(i,0);
174      ln_y(i)=log(data(i,1));
175    }
176
177    regression::Polynomial polynomialfit(2);
178    polynomialfit.fit(x,ln_y);
179    utility::vector fit=polynomialfit.fit_parameters();
180    if (fabs(fit[0]-1.012229646706 + fit[1]-0.012561322528 +
181             fit[2]+1.159674470130)>1e-11) {  // Jari, fix number!
182      *error << "regression_Polynomial: cannot reproduce fit." << std::endl;
183      ok=false;
184    }
185  }
186
187  *error << "  testing regression::PolynomialWeighted" << std::endl;
188  regression::PolynomialWeighted pol_weighted(2);
189  // ok = equal(polynomial, pol_weighted, error) && ok;
190
191  ok = multidim(error) && ok;
192
193  if (error!=&std::cerr)
194    delete error;
195
196  return (ok ? 0 : -1);
197}
198
199
200bool equal(regression::OneDimensional& r, 
201           regression::OneDimensionalWeighted& wr, 
202           std::ostream* error)
203{
204  bool ok=true;
205  utility::vector x(5); x(0)=1970; x(1)=1980; x(2)=1990; x(3)=2000; x(4)=2010;
206  utility::vector y(5); y(0)=12;   y(1)=11;   y(2)=14;   y(3)=13;   y(4)=15;
207
208  ok = unity_weights(r, wr, x, y, error) && ok;
209  ok = rescale_weights(wr, x, y, error) && ok;
210  ok = zero_weights(wr, x, y, error) && ok;
211  return ok;
212}
213
214
215bool unity_weights(regression::OneDimensional& r, 
216                   regression::OneDimensionalWeighted& rw, 
217                   const utility::vector& x, const utility::vector& y,
218                   std::ostream* error)
219{
220  *error << "    testing unity weights equal to non-weighted version.\n"; 
221  bool ok=true;
222  utility::vector w(x.size(), 1.0);
223  r.fit(x,y);
224  rw.fit(x,y,w);
225  if (r.predict(2000) != rw.predict(2000)){
226    ok = false;
227    *error << "Error: predict not equal\n" 
228           << "   weighted: " << rw.predict(2000) << "\n"
229           << "   non-weighted: " << r.predict(2000)
230           << std::endl;
231  }
232  if (fabs(r.prediction_error2(2000)-rw.prediction_error2(2000))>10e-7){
233    ok = false;
234    *error << "Error: prediction_error2 not equal non-weighted version.\n" 
235           << "       weighted: " << rw.prediction_error2(2000) << "\n"
236           << "   non-weighted: " << r.prediction_error2(2000)
237           << std::endl;
238  }
239  if (fabs(r.s2()-rw.s2())>10E-7){
240    ok = false;
241    *error << "Error: r2 not equal non-weighted version." << std::endl;
242    *error << "weighted r2 = " << rw.r2() << std::endl;
243    *error << "non-weighted r2 = " << r.r2() << std::endl;
244  }
245  if (fabs(r.s2()-rw.s2())>10E-7){
246    ok = false;
247    *error << "Error: s2 not equal non-weighted version." << std::endl;
248    *error << "weighted s2 = " << rw.s2() << std::endl;
249    *error << "non-weighted s2 = " << r.s2() << std::endl;
250  }
251  if (fabs(r.standard_error2(2000)-rw.standard_error2(2000))>10e-7){
252    ok = false;
253    *error << "Error: standard_error not equal non-weighted version." 
254           << std::endl;
255  }
256  return ok;
257} 
258
259
260bool rescale_weights(regression::OneDimensionalWeighted& wr, 
261                     const utility::vector& x, const utility::vector& y,
262                     std::ostream* error)
263{
264  *error << "    testing rescaling weights.\n"; 
265  bool ok = true;
266  utility::vector w(5);  w(0)=1.0;  w(1)=1.0;  w(2)=0.5;  w(3)=0.2;  w(4)=0.2;
267  wr.fit(x,y,w);
268  double predict = wr.predict(2000);
269  double prediction_error2 = wr.prediction_error2(2000);
270  double r2 = wr.r2();
271  double s2 = wr.s2();
272  double standard_error2 = wr.standard_error2(2000);
273
274  w.scale(2);
275  wr.fit(x,y,w);
276  if (wr.predict(2000) != predict){
277    ok = false;
278    *error << "Error: predict not equal after rescaling.\n";
279  }
280  if (wr.prediction_error2(2000,2) != prediction_error2){
281    ok = false;
282    *error << "Error: prediction_error2 not equal after rescaling.\n";
283  }
284  if (wr.r2() != r2){
285    ok = false;
286    *error << "Error: r2 not equal after rescaling.\n";
287  }
288  if (wr.s2(2) != s2){
289    ok = false;
290    *error << "Error: s2 not equal after rescaling.\n";
291    *error << "       s2 = " << s2 << " and after doubling weights.\n";
292    *error << "       s2 = " << wr.s2(2) << "\n";
293  }
294  if (wr.standard_error2(2000) != standard_error2){
295    ok = false;
296    *error << "Error: standard_error2 not equal after rescaling.\n";
297  }
298  return ok; 
299}
300
301
302bool zero_weights(regression::OneDimensionalWeighted& wr, 
303                  const utility::vector& x, const utility::vector& y,
304                  std::ostream* error)
305{
306  *error << "    testing zero weights equal to missing value.\n"; 
307  bool ok = true;
308  utility::vector w(5);  w(0)=1.0;  w(1)=1.0;  w(2)=0.5;  w(3)=0.2;  w(4)=0;
309  wr.fit(x,y,w);
310  double predict = wr.predict(2000);
311  double prediction_error2 = wr.prediction_error2(2000);
312  double r2 = wr.r2();
313  double s2 = wr.s2();
314  double standard_error2 = wr.standard_error2(2000);
315
316  utility::vector x2(4);
317  utility::vector y2(4);
318  utility::vector w2(4);
319  for (size_t i=0; i<4; ++i){
320    x2(i) = x(i);
321    y2(i) = y(i);
322    w2(i) = w(i);
323  }
324
325  wr.fit(x2,y2,w2);
326  if (wr.predict(2000) != predict){
327    ok = false;
328    *error << "Error: predict not equal.\n";
329  }
330  if (wr.prediction_error2(2000) != prediction_error2){
331    ok = false;
332    *error << "Error: prediction_error2 not equal.\n";
333  }
334  if (wr.r2() != r2){
335    ok = false;
336    *error << "Error: r2 not equal.\n";
337  }
338  if (wr.s2() != s2){
339    ok = false;
340    *error << "Error: s2 not equal.\n";
341  }
342  if (wr.standard_error2(2000) != standard_error2){
343    ok = false;
344    *error << "Error: standard_error2 not equal.\n";
345  }
346  return ok; 
347}
348
349
350bool multidim(std::ostream* error)
351{
352  bool ok = true;
353  *error << "  testing regression::MultiDimensionalWeighted" << std::endl;
354  utility::vector x(5); x(0)=1970; x(1)=1980; x(2)=1990; x(3)=2000; x(4)=2010;
355  utility::vector y(5); y(0)=12;   y(1)=11;   y(2)=14;   y(3)=13;   y(4)=15;
356  utility::vector w(5,1.0);
357 
358  utility::matrix data(5,3);
359  for (size_t i=0; i<data.rows(); ++i){
360    data(i,0)=1;
361    data(i,1)=x(i);
362    data(i,2)=x(i)*x(i);
363  }
364  regression::MultiDimensional md;
365  md.fit(data,y);
366  regression::MultiDimensionalWeighted mdw;
367  mdw.fit(data,y,w);
368  utility::vector z(3,1);
369  z(1)=2000;
370  z(2)=2000*2000;
371  if (md.predict(z) != mdw.predict(z)){
372    ok = false;
373    *error << "Error: predict not equal\n" 
374           << "   weighted: " << mdw.predict(1) << "\n"
375           << "   non-weighted: " << md.predict(1)
376           << std::endl;
377  }
378  /*
379  if (md.standard_error2(z) != mdw.standard_error2(z)){
380    ok = false;
381    *error << "Error: standard_error2 not equal\n"
382           << "   weighted: " << mdw.standard_error2(z) << "\n"
383           << "   non-weighted: " << md.standard_error2(z)
384           << std::endl;
385  }
386  */
387  return ok;
388}
389
390
391bool Local_test(regression::OneDimensionalWeighted& r, 
392                regression::Kernel& k)
393{
394  regression::Local rl(r,k);
395  for (size_t i=0; i<1000; i++){
396    rl.add(i, 10);
397  }
398
399  rl.fit(10, 100);
400
401  utility::vector y = rl.y_predicted();
402  for (size_t i=0; i<y.size(); i++) 
403    if (y(i)!=10.0){
404      return false;
405    }
406  return true;
407}
Note: See TracBrowser for help on using the repository browser.