test/regression.cc

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

Generated on Sat Oct 14 02:21:26 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)