source: trunk/test/regression_test.cc @ 729

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

Fixes #159. Also removed some inlines in OneDimensionalWeighted? by adding source file. Refs #81

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 10.3 KB
Line 
1// $Id: regression_test.cc 729 2007-01-05 16:00:15Z 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  linear_w.fit(x,y,w);
123  double y_predicted = linear_w.predict(1990);
124  double y_predicted_err = linear_w.prediction_error2(1990);
125  if (fabs(y_predicted-12.8)>0.001){
126    *error << "error: cannot reproduce fit." << std::endl;
127    *error << "predicted value: " << y_predicted << " expected 12.8" 
128           << std::endl;
129    ok=false;
130  }
131
132  // testing regression::NaiveWeighted
133  *error << "  testing regression::NaiveWeighted" << std::endl;
134  regression::NaiveWeighted naive_w;
135  regression::Naive naive;
136  ok = ok && equal(naive, naive_w, error);
137  naive_w.fit(x,y,w);
138
139  y_predicted=naive_w.predict(0.0);
140  y_predicted_err=naive_w.prediction_error2(0.0);
141  if (y_predicted!=(0.1*12+0.2*11+0.3*14+0.4*13)) {
142    *error << "regression_NaiveWeighted: cannot reproduce fit." << std::endl;
143    *error << "returned value: " << y_predicted << std::endl;
144    *error << "expected: " << 0.1*12+0.2*11+0.3*14+0.4*13 << std::endl;
145    ok=false;
146  }
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::Linear" << std::endl;
186  regression::Linear lin;
187 
188  *error << "  testing regression::PolynomialWeighted" << std::endl;
189  regression::PolynomialWeighted pol_weighted(2);
190
191  if (error!=&std::cerr)
192    delete error;
193
194  return (ok ? 0 : -1);
195}
196
197
198bool equal(regression::OneDimensional& r, 
199           regression::OneDimensionalWeighted& wr, 
200           std::ostream* error)
201{
202  bool ok=true;
203  utility::vector x(5); x(0)=1970; x(1)=1980; x(2)=1990; x(3)=2000; x(4)=2010;
204  utility::vector y(5); y(0)=12;   y(1)=11;   y(2)=14;   y(3)=13;   y(4)=15;
205
206  ok = unity_weights(r, wr, x, y, error) && ok;
207  ok = rescale_weights(wr, x, y, error) && ok;
208  ok = zero_weights(wr, x, y, error) && ok;
209  return ok;
210}
211
212
213bool unity_weights(regression::OneDimensional& r, 
214                   regression::OneDimensionalWeighted& rw, 
215                   const utility::vector& x, const utility::vector& y,
216                   std::ostream* error)
217{
218  *error << "    testing unity weights equal to non-weighted version.\n"; 
219  bool ok=true;
220  utility::vector w(x.size(), 1.0);
221  r.fit(x,y);
222  rw.fit(x,y,w);
223  if (r.predict(2000) != rw.predict(2000)){
224    ok = false;
225    *error << "Error: predict not equal" << std::endl;
226  }
227  if (r.prediction_error2(2000) != rw.prediction_error2(2000)){
228    ok = false;
229    *error << "Error: prediction_error2 not equal non-weighted version." 
230           << std::endl;
231  }
232  if (r.r2() != rw.r2()){
233    ok = false;
234    *error << "Error: r2 not equal non-weighted version." << std::endl;
235  }
236  if (r.s2() != rw.s2()){
237    ok = false;
238    *error << "Error: s2 not equal non-weighted version." << std::endl;
239  }
240  if (r.standard_error2(2000) != rw.standard_error2(2000)){
241    ok = false;
242    *error << "Error: standard_error not equal non-weighted version." 
243           << std::endl;
244  }
245  return ok;
246} 
247
248
249bool rescale_weights(regression::OneDimensionalWeighted& wr, 
250                     const utility::vector& x, const utility::vector& y,
251                     std::ostream* error)
252{
253  *error << "    testing rescaling weights.\n"; 
254  bool ok = true;
255  utility::vector w(5);  w(0)=1.0;  w(1)=1.0;  w(2)=0.5;  w(3)=0.2;  w(4)=0.2;
256  wr.fit(x,y,w);
257  double predict = wr.predict(2000);
258  double prediction_error2 = wr.prediction_error2(2000);
259  double r2 = wr.r2();
260  double s2 = wr.s2();
261  double standard_error2 = wr.standard_error2(2000);
262
263  w.scale(2);
264  wr.fit(x,y,w);
265  if (wr.predict(2000) != predict){
266    ok = false;
267    *error << "Error: predict not equal after rescaling.\n";
268  }
269  if (wr.prediction_error2(2000,2) != prediction_error2){
270    ok = false;
271    *error << "Error: prediction_error2 not equal after rescaling.\n";
272  }
273  if (wr.r2() != r2){
274    ok = false;
275    *error << "Error: r2 not equal after rescaling.\n";
276  }
277  if (wr.s2(2) != s2){
278    ok = false;
279    *error << "Error: s2 not equal after rescaling.\n";
280    *error << "       s2 = " << s2 << " and after doubling weights.\n";
281    *error << "       s2 = " << wr.s2() << "\n";
282  }
283  if (wr.standard_error2(2000) != standard_error2){
284    ok = false;
285    *error << "Error: standard_error2 not equal after rescaling.\n";
286  }
287  return ok; 
288}
289
290
291bool zero_weights(regression::OneDimensionalWeighted& wr, 
292                  const utility::vector& x, const utility::vector& y,
293                  std::ostream* error)
294{
295  *error << "    testing zero weights equal to missing value.\n"; 
296  bool ok = true;
297  utility::vector w(5);  w(0)=1.0;  w(1)=1.0;  w(2)=0.5;  w(3)=0.2;  w(4)=0;
298  wr.fit(x,y,w);
299  double predict = wr.predict(2000);
300  double prediction_error2 = wr.prediction_error2(2000);
301  double r2 = wr.r2();
302  double s2 = wr.s2();
303  double standard_error2 = wr.standard_error2(2000);
304
305  utility::vector x2(4);
306  utility::vector y2(4);
307  utility::vector w2(4);
308  for (size_t i=0; i<4; ++i){
309    x2(i) = x(i);
310    y2(i) = y(i);
311    w2(i) = w(i);
312  }
313
314  wr.fit(x2,y2,w2);
315  if (wr.predict(2000) != predict){
316    ok = false;
317    *error << "Error: predict not equal.\n";
318  }
319  if (wr.prediction_error2(2000) != prediction_error2){
320    ok = false;
321    *error << "Error: prediction_error2 not equal.\n";
322  }
323  if (wr.r2() != r2){
324    ok = false;
325    *error << "Error: r2 not equal.\n";
326  }
327  if (wr.s2() != s2){
328    ok = false;
329    *error << "Error: s2 not equal.\n";
330  }
331  if (wr.standard_error2(2000) != standard_error2){
332    ok = false;
333    *error << "Error: standard_error2 not equal.\n";
334  }
335  return ok; 
336}
337
338
339bool Local_test(regression::OneDimensionalWeighted& r, 
340                regression::Kernel& k)
341{
342  regression::Local rl(r,k);
343  for (size_t i=0; i<1000; i++){
344    rl.add(i, 10);
345  }
346
347  rl.fit(10, 100);
348
349  utility::vector y = rl.y_predicted();
350  for (size_t i=0; i<y.size(); i++) 
351    if (y(i)!=10.0){
352      return false;
353    }
354  return true;
355}
Note: See TracBrowser for help on using the repository browser.