source: trunk/test/regression.cc @ 3114

Last change on this file since 3114 was 3114, checked in by peter, 4 years ago

update copyright years

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 17.2 KB
Line 
1// $Id: regression.cc 3114 2013-11-10 23:51:47Z peter $
2
3/*
4        Copyright (C) 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
5        Copyright (C) 2009 Peter Johansson
6        Copyright (C) 2012 Jari Häkkinen, Peter Johansson
7        Copyright (C) 2013 Peter Johansson
8
9        This file is part of the yat library, http://dev.thep.lu.se/yat
10
11        The yat library is free software; you can redistribute it and/or
12        modify it under the terms of the GNU General Public License as
13        published by the Free Software Foundation; either version 3 of the
14        License, or (at your option) any later version.
15
16        The yat library is distributed in the hope that it will be useful,
17        but WITHOUT ANY WARRANTY; without even the implied warranty of
18        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19        General Public License for more details.
20
21        You should have received a copy of the GNU General Public License
22        along with yat. If not, see <http://www.gnu.org/licenses/>.
23*/
24
25#include <config.h>
26
27#include "Suite.h"
28
29#include "yat/regression/KernelBox.h"
30#include "yat/regression/Linear.h"
31#include "yat/regression/LinearWeighted.h"
32#include "yat/regression/Local.h"
33#include "yat/regression/Naive.h"
34#include "yat/regression/NaiveWeighted.h"
35#include "yat/regression/Polynomial.h"
36#include "yat/regression/PolynomialWeighted.h"
37#include "yat/utility/Matrix.h"
38#include "yat/utility/Vector.h"
39
40#include <cmath>
41
42#include <fstream>
43#include <iostream>
44
45
46using namespace theplu::yat;
47
48void equal(regression::OneDimensional&, regression::OneDimensionalWeighted&, 
49                                         test::Suite&); 
50
51void multidim(test::Suite& suite);
52
53void unity_weights(regression::OneDimensional&, 
54                                                                         regression::OneDimensionalWeighted&, 
55                                                                         const utility::Vector&, const utility::Vector&,
56                                                                         test::Suite&); 
57
58void rescale_weights(regression::OneDimensionalWeighted&, 
59                                                                                 const utility::Vector&, const utility::Vector&,
60                                                                                 test::Suite&); 
61
62void zero_weights(regression::OneDimensionalWeighted&, 
63                                                                        const utility::Vector&, const utility::Vector&,
64                                                                        test::Suite&); 
65
66
67bool Local_test(regression::OneDimensionalWeighted&, 
68                                                                regression::Kernel&, test::Suite&);
69
70int main(int argc, char* argv[])
71{
72        test::Suite suite(argc, argv);
73
74        suite.err() << "  testing regression" << std::endl;
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 ( !suite.equal(linear.beta(),polynomial.fit_parameters()(1), 1000) ){
87                suite.err() << "error: beta and fit_parameters(1) not equal" << std::endl;
88                suite.err() << "       beta = " << linear.beta() << std::endl;
89                suite.err() << "       fit_parameters(1) = " 
90                                         << polynomial.fit_parameters()(1) << std::endl;
91                suite.add(false);
92        }
93        if (!suite.equal(polynomial.fit_parameters()(0),
94                                                                         linear.alpha()-linear.beta()*1985, 10000)){
95                suite.err() << "error: fit_parameters(0) = " 
96                                         << polynomial.fit_parameters()(0)<< std::endl;
97                suite.err() << "error: alpha-beta*m_x = " 
98                                         << linear.alpha()-linear.beta()*1985 << std::endl;
99                suite.add(false);
100        }
101        if ( !suite.equal(polynomial.chisq(), linear.chisq(), 100) ){
102                suite.err() << "error: chisq not same in linear and polynomial(1)" 
103                                         << std::endl;
104                suite.add(false);
105        }
106        if ( !suite.equal(polynomial.predict(1.0),linear.predict(1.0), 1000) ){
107                suite.err() << "error: predict not same in linear and polynomial(1)" 
108                                         << std::endl;
109                suite.add(false);
110        }
111        if (!suite.equal(polynomial.standard_error2(1985),
112                                                                         linear.standard_error2(1985), 100000)){
113                suite.err() << "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                suite.add(false);
118        }
119
120        suite.err() << "  testing regression::LinearWeighted" << std::endl;
121        regression::LinearWeighted linear_w;
122        equal(linear, linear_w, suite);
123        linear_w.fit(x,y,w);
124        double y_predicted = linear_w.predict(1990);
125        if (!suite.equal(y_predicted,12.8) ){
126                suite.err() << "error: cannot reproduce fit." << std::endl;
127                suite.err() << "predicted value: " << y_predicted << " expected 12.8" 
128                                         << std::endl;
129                suite.add(false);
130        }
131        double y_predicted_err = linear_w.prediction_error2(1990);
132        if (!suite.equal_sqrt(y_predicted_err,1.2) ){
133                suite.err() << "error: cannot reproduce fit." << std::endl;
134                suite.err() << "predicted error: " << y_predicted_err << " expected 1.2\n";
135                suite.add(false);
136        }
137
138        // Comparing LinearWeighted and PolynomialWeighted(1)
139        suite.err() << "    comparing LinearWeighted and PolynomialWeighted(1)"
140                                                        << std::endl;
141        linear_w.fit(x,y,w);
142        regression::PolynomialWeighted polynomial_w(1);
143        polynomial_w.fit(x,y,w);
144        if ( !suite.equal(linear_w.beta(), polynomial_w.fit_parameters()(1),10000) ){
145                suite.err() << "error: beta and fit_parameters(1) not equal" << std::endl;
146                suite.err() << "       beta = " << linear_w.beta() << std::endl;
147                suite.err() << "       fit_parameters(1) = " 
148                                         << polynomial_w.fit_parameters()(1) << std::endl;
149                suite.add(false);
150        }
151        if ( !suite.equal(polynomial_w.fit_parameters()(0),
152                                                                                linear_w.alpha()-linear_w.beta()*1990, 10000) ){
153                suite.err() << "error: fit_parameters(0) = " 
154                                         << polynomial.fit_parameters()(0)<< std::endl;
155                suite.err() << "error: alpha-beta*m_x = " 
156                                         << linear.alpha()-linear.beta()*1990 << std::endl;
157                suite.add(false);
158        }
159        if ( !suite.equal(polynomial_w.s2(),linear_w.s2(), 10) ){
160                suite.err() << "error: chisq not same in linear and polynomial(1)" 
161                                         << std::endl;
162                suite.add(false);
163        }
164        if ( !suite.equal(polynomial_w.predict(1.0), linear_w.predict(1.0), 10000) ){
165                suite.err() << "error: predict not same in linear and polynomial(1)" 
166                                         << std::endl;
167                suite.add(false);
168        }
169        if ( !suite.equal(polynomial_w.standard_error2(1985), 
170                                                                                linear_w.standard_error2(1985), 100000) ){
171                suite.err() << "error: standard_error not same in linear and polynomial(1)" 
172                                         << "\n  polynomial: " << polynomial_w.standard_error2(1985)
173                                         << "\n  linear: " << linear_w.standard_error2(1985)
174                                         << std::endl;
175                suite.add(false);
176        }
177
178        // testing regression::NaiveWeighted
179        suite.err() << "  testing regression::NaiveWeighted" << std::endl;
180        regression::NaiveWeighted naive_w;
181        regression::Naive naive;
182        equal(naive, naive_w, suite);
183        naive_w.fit(x,y,w);
184
185        y_predicted=naive_w.predict(0.0);
186        y_predicted_err=naive_w.prediction_error2(0.0);
187        if (!suite.equal(y_predicted,0.1*12+0.2*11+0.3*14+0.4*13)) {
188                suite.err() << "regression_NaiveWeighted: cannot reproduce fit.\n";
189                suite.err() << "returned value: " << y_predicted << std::endl;
190                suite.err() << "expected: " << 0.1*12+0.2*11+0.3*14+0.4*13 << std::endl;
191                suite.add(false);
192        }
193
194        // testing regression::Local
195        suite.err() << "  testing regression::Local" << std::endl;
196        regression::KernelBox kb;
197        regression::LinearWeighted rl;
198        if (!Local_test(rl,kb, suite)) {
199                suite.err() << "regression_Local: Linear cannot reproduce fit." << std::endl;
200                suite.add(false);
201        }
202        regression::NaiveWeighted rn;
203        if (!Local_test(rn,kb, suite)) {
204                suite.err() << "regression_Local: Naive cannot reproduce fit." << std::endl;
205                suite.add(false);
206        }
207
208        // testing regression::Polynomial
209        suite.err() << "  testing regression::Polynomial" << std::endl;
210        {
211                std::ifstream s(test::filename("data/regression_gauss.data").c_str());
212                utility::Matrix data(s);
213                utility::Vector x(data.rows());
214                utility::Vector ln_y(data.rows());
215                for (size_t i=0; i<data.rows(); ++i) {
216                        x(i)=data(i,0);
217                        ln_y(i)=log(data(i,1));
218                }
219
220                regression::Polynomial polynomialfit(2);
221                polynomialfit.fit(x,ln_y);
222                utility::Vector fit=polynomialfit.fit_parameters();
223                suite.add(suite.equal_fix(fit(0), 1.012229646706, 1e-11));
224                suite.add(suite.equal_fix(fit(1), 0.012561322528, 1e-11));
225                suite.add(suite.equal_fix(fit(2), -1.159674470130, 1e-11));
226        }
227
228        suite.err() << "  testing regression::PolynomialWeighted" << std::endl;
229        regression::PolynomialWeighted pol_weighted(2);
230        regression::Polynomial polynomial2(2);
231        equal(polynomial2, pol_weighted, suite);
232
233        multidim(suite);
234
235        return suite.return_value();
236}
237
238
239void equal(regression::OneDimensional& r, 
240                                         regression::OneDimensionalWeighted& wr, 
241                                         test::Suite& suite)
242{
243        utility::Vector x(5); x(0)=1970; x(1)=1980; x(2)=1990; x(3)=2000; x(4)=2010;
244        utility::Vector y(5); y(0)=12;   y(1)=11;   y(2)=14;   y(3)=13;   y(4)=15;
245
246        unity_weights(r, wr, x, y, suite);
247        rescale_weights(wr, x, y, suite);
248        zero_weights(wr, x, y, suite);
249}
250
251
252void unity_weights(regression::OneDimensional& r, 
253                                                                         regression::OneDimensionalWeighted& rw, 
254                                                                         const utility::Vector& x, const utility::Vector& y,
255                                                                         test::Suite& suite)
256{
257        suite.err() << "    testing unity weights equal to non-weighted version.\n"; 
258        utility::Vector w(x.size(), 1.0);
259        r.fit(x,y);
260        rw.fit(x,y,w);
261        if (!suite.equal(r.predict(2000), rw.predict(2000)) ) {
262                suite.add(false);
263                suite.err() << "Error: predict not equal\n" 
264                                         << "   weighted: " << rw.predict(2000) << "\n"
265                                         << "   non-weighted: " << r.predict(2000)
266                                         << std::endl;
267        }
268        if (!suite.equal(r.s2(), rw.s2(1.0)) ){
269                suite.add(false);
270                suite.err() << "Error: s2 not equal non-weighted version." << std::endl;
271                suite.err() << "weighted s2 = " << rw.s2(1.0) << std::endl;
272                suite.err() << "non-weighted s2 = " << r.s2() << std::endl;
273        }
274        if (!suite.equal_sqrt(r.standard_error2(2000), rw.standard_error2(2000), 20)){
275                suite.add(false);
276                suite.err() << "Error: standard_error not equal non-weighted version." 
277                                         << std::endl;
278        }
279        if (!suite.equal(r.r2(), rw.r2()) ){
280                suite.add(false);
281                suite.err() << "Error: r2 not equal non-weighted version." << std::endl;
282                suite.err() << "weighted r2 = " << rw.r2() << std::endl;
283                suite.err() << "non-weighted r2 = " << r.r2() << std::endl;
284        }
285        if (!suite.equal_sqrt(r.prediction_error2(2000), rw.prediction_error2(2000,1),
286                                                                                                100) ){
287                suite.add(false);
288                suite.err() << "Error: prediction_error2 not equal non-weighted version.\n" 
289                                         << "       weighted: " << rw.prediction_error2(2000,1) << "\n"
290                                         << "   non-weighted: " << r.prediction_error2(2000)
291                                         << std::endl;
292        }
293} 
294
295
296void rescale_weights(regression::OneDimensionalWeighted& wr, 
297                                                                                 const utility::Vector& x, const utility::Vector& y,
298                                                                                 test::Suite& suite)
299{
300        suite.err() << "    testing rescaling weights.\n"; 
301        utility::Vector w(5);  w(0)=1.0;  w(1)=1.0;  w(2)=0.5;  w(3)=0.2;  w(4)=0.2;
302        wr.fit(x,y,w);
303        double predict = wr.predict(2000);
304        double prediction_error2 = wr.prediction_error2(2000);
305        double r2 = wr.r2();
306        double s2 = wr.s2();
307        double standard_error2 = wr.standard_error2(2000);
308
309        w*=2;
310        wr.fit(x,y,w);
311        if (!suite.equal(wr.predict(2000), predict, 10000) ){
312                suite.add(false);
313                suite.err() << "Error: predict not equal after rescaling.\n";
314                suite.err() << "       predict = " << predict
315                                         << " and after doubling weights.\n";
316                suite.err() << "       predict = " << wr.predict(2000) << "\n";
317        }
318        if (!suite.equal(wr.s2(2), s2, 14000) ){
319                suite.add(false);
320                suite.err() << "Error: s2 not equal after rescaling.\n";
321                suite.err() << "       s2 = " << s2 << " and after doubling weights.\n";
322                suite.err() << "       s2 = " << wr.s2(2) << "\n";
323                suite.err() << "difference " << s2-wr.s2(2.0) << std::endl;
324        }
325        if (!suite.equal_sqrt(wr.standard_error2(2000), standard_error2, 100) ){
326                suite.add(false);
327                suite.err() << "Error: standard_error2 not equal after rescaling.\n";
328                suite.err() << "       standard_error2 = " << standard_error2
329                                         << " and after doubling weights.\n";
330                suite.err() << "       standard_error2 = " 
331                                         << wr.standard_error2(2000) << "\n";
332                suite.err() << " difference: " << wr.standard_error2(2000)-standard_error2
333                                         << std::endl;
334        }
335        if (!suite.equal(wr.r2(), r2, 10000) ){
336                suite.add(false);
337                suite.err() << "Error: r2 not equal after rescaling.\n";
338        }
339        if (!suite.equal_sqrt(wr.prediction_error2(2000,2), prediction_error2, 20) ){
340                suite.add(false);
341                suite.err() << "Error: prediction_error2 not equal after rescaling.\n";
342                suite.err() << "       prediction_error2 = " << prediction_error2
343                                         << " and after doubling weights.\n";
344                suite.err() << "       prediction_error2 = " 
345                                         << wr.prediction_error2(2000,2) << "\n";
346        }
347}
348
349
350void zero_weights(regression::OneDimensionalWeighted& wr, 
351                                                                        const utility::Vector& x, const utility::Vector& y,
352                                                                        test::Suite& suite)
353{
354        suite.err() << "    testing zero weights equal to missing value.\n"; 
355        utility::Vector w(5);  w(0)=1.0;  w(1)=1.0;  w(2)=0.5;  w(3)=0.2;  w(4)=0;
356        wr.fit(x,y,w);
357        double predict = wr.predict(2000);
358        double prediction_error2 = wr.prediction_error2(2000);
359        double r2 = wr.r2();
360        double s2 = wr.s2();
361        double standard_error2 = wr.standard_error2(2000);
362
363        utility::Vector x2(4);
364        utility::Vector y2(4);
365        utility::Vector w2(4);
366        for (size_t i=0; i<4; ++i){
367                x2(i) = x(i);
368                y2(i) = y(i);
369                w2(i) = w(i);
370        }
371
372        wr.fit(x2,y2,w2);
373        if (!suite.equal(wr.predict(2000), predict, 10000) ) {
374                suite.add(false);
375                suite.err() << "Error: predict not equal.\n";
376                suite.err() << "       weighted predict: " << wr.predict(2000) << "\n";
377                suite.err() << "       unweighted predict: " << predict << "\n";
378                suite.err() << "       difference: " << wr.predict(2000)-predict << "\n";
379
380        }
381        if (!suite.equal(wr.prediction_error2(2000), prediction_error2) ) {
382                suite.add(false);
383                suite.err() << "Error: prediction_error2 not equal.\n";
384        }
385        if (!suite.equal(wr.r2(), r2) ) {
386                suite.add(false);
387                suite.err() << "Error: r2 not equal.\n";
388                suite.err() << "   r2: " << r2 << "\n";
389                suite.err() << "   r2: " << wr.r2() << "\n";
390        }
391        if (!suite.equal(wr.s2(), s2) ) {
392                suite.add(false);
393                suite.err() << "Error: s2 not equal.\n";
394        }
395        if (!suite.equal(wr.standard_error2(2000), standard_error2) ) {
396                suite.add(false);
397                suite.err() << "Error: standard_error2 not equal.\n";
398        }
399}
400
401
402void multidim(test::Suite& suite)
403{
404        suite.err() << "  testing regression::MultiDimensionalWeighted" << std::endl;
405        utility::Vector x(5); x(0)=1970; x(1)=1980; x(2)=1990; x(3)=2000; x(4)=2010;
406        utility::Vector y(5); y(0)=12;   y(1)=11;   y(2)=14;   y(3)=13;   y(4)=15;
407        utility::Vector w(5,1.0);
408       
409        utility::Matrix data(5,3);
410        for (size_t i=0; i<data.rows(); ++i){
411                data(i,0)=1;
412                data(i,1)=x(i);
413                data(i,2)=x(i)*x(i);
414        }
415        regression::MultiDimensional md;
416        md.fit(data,y);
417        regression::MultiDimensionalWeighted mdw;
418        mdw.fit(data,y,w);
419        utility::Vector z(3,1);
420        z(1)=2000;
421        z(2)=2000*2000;
422        if (!suite.equal(md.predict(z), mdw.predict(z))){
423                suite.add(false);
424                suite.err() << "Error: predict not equal\n" 
425                                         << "   weighted: " << mdw.predict(z) << "\n"
426                                         << "   non-weighted: " << md.predict(z)
427                                         << std::endl;
428        }
429
430        if (!suite.equal_sqrt(md.standard_error2(z), mdw.standard_error2(z),20) ){
431                suite.add(false);
432                suite.err() << "Error: standard_error2 not equal\n" 
433                                         << "   weighted: " << mdw.standard_error2(z) << "\n"
434                                         << "   non-weighted: " << md.standard_error2(z)
435                                         << std::endl;
436        }
437        if (!suite.equal_sqrt(md.prediction_error2(z), mdw.prediction_error2(z,1.0),
438                                                                                                20) ){
439                suite.add(false);
440                suite.err() << "Error: prediction_error2 not equal\n" 
441                                         << "   weighted: " << mdw.prediction_error2(z,1.0) << "\n"
442                                         << "   non-weighted: " << md.prediction_error2(z)
443                                         << std::endl;
444        }
445
446        w(0)=1.0;  w(1)=1.0;  w(2)=0.5;  w(3)=0.2;  w(4)=0.2;
447        mdw.fit(data,y,w);
448        double predict = mdw.predict(z);
449        double prediction_error2 = mdw.prediction_error2(z, 1.0);
450        double s2 = mdw.s2(1.0);
451        double standard_error2 = mdw.standard_error2(z);
452
453        w*=2;
454        mdw.fit(data,y,w);
455        if (!suite.equal(mdw.predict(z), predict, 10000) ){
456                suite.add(false);
457                suite.err() << "Error: predict not equal after rescaling.\n";
458                suite.err() << "   predict = " << predict
459                                                                << " and after doubling weights.\n";
460                suite.err() << "   predict = " << mdw.predict(z) << "\n";
461        }
462        if (!suite.equal_sqrt(mdw.prediction_error2(z,2), prediction_error2,20) ){ 
463                suite.add(false);
464                suite.err() << "Error: prediction_error2 not equal after rescaling.\n";
465                suite.err() << "       predict_error2 = " << prediction_error2
466                                         << " and after doubling weights.\n";
467                suite.err() << "       predict_error2 = " 
468                                                                << mdw.prediction_error2(z,2) << "\n";
469        }
470        if (!suite.equal(mdw.s2(2), s2, 14000) ){
471                suite.add(false);
472                suite.err() << "Error: s2 not equal after rescaling.\n";
473                suite.err() << "       s2 = " << s2 << " and after doubling weights.\n";
474                suite.err() << "       s2 = " << mdw.s2(2) << "\n";
475        }
476        if (!suite.equal_sqrt(mdw.standard_error2(z), standard_error2, 100) ){
477                suite.add(false);
478                suite.err() << "Error: standard_error2 not equal after rescaling.\n";
479                suite.err() << " standard_error2 = " << standard_error2
480                                         << " and after doubling weights.\n";
481                suite.err() << " standard_error2 = " << mdw.standard_error2(z) << "\n";
482        }
483}
484
485
486bool Local_test(regression::OneDimensionalWeighted& r, 
487                                                                regression::Kernel& k, test::Suite& suite)
488{
489        regression::Local rl(r,k);
490        for (size_t i=0; i<1000; i++){
491                rl.add(i, 10);
492        }
493
494        rl.fit(10, 100);
495        if (rl.x().size()!=1000/10)
496                return false;
497        for (size_t i=0; i+1<rl.x().size(); ++i)
498                if (!suite.equal(rl.x()(i+1)-rl.x()(i),10.0))
499                        return false;
500
501        utility::Vector y(rl.y_predicted());
502        for (size_t i=0; i<y.size(); i++) 
503                if (!suite.equal(y(i),10.0)){
504                        return false;
505                }
506
507        return true;
508}
Note: See TracBrowser for help on using the repository browser.