SHOGUN
v3.2.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 2013 Roman Votyakov 00008 * Copyright (C) 2012 Jacob Walker 00009 * Copyright (C) 2013 Roman Votyakov 00010 * 00011 * Adapted from the GPML toolbox, specifically likT.m 00012 * http://www.gaussianprocess.org/gpml/code/matlab/doc/ 00013 */ 00014 00015 #include <shogun/machine/gp/StudentsTLikelihood.h> 00016 00017 #ifdef HAVE_EIGEN3 00018 00019 #include <shogun/mathematics/Function.h> 00020 #include <shogun/mathematics/Integration.h> 00021 #include <shogun/labels/RegressionLabels.h> 00022 #include <shogun/mathematics/Statistics.h> 00023 #include <shogun/mathematics/Math.h> 00024 #include <shogun/mathematics/eigen3.h> 00025 00026 using namespace shogun; 00027 using namespace Eigen; 00028 00029 namespace shogun 00030 { 00031 00032 #ifndef DOXYGEN_SHOULD_SKIP_THIS 00033 00035 class CNormalPDF : public CFunction 00036 { 00037 public: 00039 CNormalPDF() 00040 { 00041 m_mu=0.0; 00042 m_sigma=1.0; 00043 } 00044 00050 CNormalPDF(float64_t mu, float64_t sigma) 00051 { 00052 m_mu=mu; 00053 m_sigma=sigma; 00054 } 00055 00060 void set_mu(float64_t mu) { m_mu=mu; } 00061 00066 void set_sigma(float64_t sigma) { m_sigma=sigma; } 00067 00074 virtual float64_t operator() (float64_t x) 00075 { 00076 return (1.0/(CMath::sqrt(2*CMath::PI)*m_sigma))* 00077 CMath::exp(-CMath::sq(x-m_mu)/(2.0*CMath::sq(m_sigma))); 00078 } 00079 00080 private: 00081 /* mean */ 00082 float64_t m_mu; 00083 00084 /* standard deviation */ 00085 float64_t m_sigma; 00086 }; 00087 00091 class CStudentsTPDF : public CFunction 00092 { 00093 public: 00095 CStudentsTPDF() 00096 { 00097 m_sigma=1.0; 00098 m_mu=0.0; 00099 m_nu=3.0; 00100 } 00101 00108 CStudentsTPDF(float64_t sigma, float64_t nu, float64_t mu) 00109 { 00110 m_sigma=sigma; 00111 m_mu=mu; 00112 m_nu=nu; 00113 } 00114 00119 void set_sigma(float64_t sigma) { m_sigma=sigma; } 00120 00125 void set_mu(float64_t mu) { m_mu=mu; } 00126 00131 void set_nu(float64_t nu) { m_nu=nu; } 00132 00140 virtual float64_t operator() (float64_t x) 00141 { 00142 float64_t lZ=CStatistics::lgamma((m_nu+1.0)/2.0)-CStatistics::lgamma(m_nu/2.0)- 00143 CMath::log(m_nu*CMath::PI*CMath::sq(m_sigma))/2.0; 00144 return CMath::exp(lZ-((m_nu+1.0)/2.0)*CMath::log(1.0+CMath::sq(x-m_mu)/ 00145 (m_nu*CMath::sq(m_sigma)))); 00146 } 00147 00148 private: 00150 float64_t m_sigma; 00151 00153 float64_t m_mu; 00154 00156 float64_t m_nu; 00157 }; 00158 00160 class CProductFunction : public CFunction 00161 { 00162 public: 00168 CProductFunction(CFunction* f, CFunction* g) 00169 { 00170 SG_REF(f); 00171 SG_REF(g); 00172 m_f=f; 00173 m_g=g; 00174 } 00175 00176 virtual ~CProductFunction() 00177 { 00178 SG_UNREF(m_f); 00179 SG_UNREF(m_g); 00180 } 00181 00188 virtual float64_t operator() (float64_t x) 00189 { 00190 return (*m_f)(x)*(*m_g)(x); 00191 } 00192 00193 private: 00195 CFunction* m_f; 00197 CFunction* m_g; 00198 }; 00199 00201 class CLinearFunction : public CFunction 00202 { 00203 public: 00205 CLinearFunction() { } 00206 00207 virtual ~CLinearFunction() { } 00208 00215 virtual float64_t operator() (float64_t x) 00216 { 00217 return x; 00218 } 00219 }; 00220 00222 class CQuadraticFunction : public CFunction 00223 { 00224 public: 00226 CQuadraticFunction() { } 00227 00228 virtual ~CQuadraticFunction() { } 00229 00236 virtual float64_t operator() (float64_t x) 00237 { 00238 return CMath::sq(x); 00239 } 00240 }; 00241 00242 #endif /* DOXYGEN_SHOULD_SKIP_THIS */ 00243 00244 CStudentsTLikelihood::CStudentsTLikelihood() : CLikelihoodModel() 00245 { 00246 init(); 00247 } 00248 00249 CStudentsTLikelihood::CStudentsTLikelihood(float64_t sigma, float64_t df) 00250 : CLikelihoodModel() 00251 { 00252 REQUIRE(sigma>0.0, "Scale parameter must be greater than zero\n") 00253 REQUIRE(df>1.0, "Number of degrees of freedom must be greater than one\n") 00254 00255 init(); 00256 m_sigma=sigma; 00257 m_df=df; 00258 } 00259 00260 void CStudentsTLikelihood::init() 00261 { 00262 m_sigma=1.0; 00263 m_df=3.0; 00264 SG_ADD(&m_df, "df", "Degrees of freedom", MS_AVAILABLE, GRADIENT_AVAILABLE); 00265 SG_ADD(&m_sigma, "sigma", "Scale parameter", MS_AVAILABLE, GRADIENT_AVAILABLE); 00266 } 00267 00268 CStudentsTLikelihood::~CStudentsTLikelihood() 00269 { 00270 } 00271 00272 CStudentsTLikelihood* CStudentsTLikelihood::obtain_from_generic( 00273 CLikelihoodModel* lik) 00274 { 00275 ASSERT(lik!=NULL); 00276 00277 if (lik->get_model_type()!=LT_STUDENTST) 00278 SG_SERROR("Provided likelihood is not of type CStudentsTLikelihood!\n") 00279 00280 SG_REF(lik); 00281 return (CStudentsTLikelihood*)lik; 00282 } 00283 00284 SGVector<float64_t> CStudentsTLikelihood::get_predictive_means( 00285 SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const 00286 { 00287 return SGVector<float64_t>(mu); 00288 } 00289 00290 SGVector<float64_t> CStudentsTLikelihood::get_predictive_variances( 00291 SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const 00292 { 00293 SGVector<float64_t> result(s2); 00294 Map<VectorXd> eigen_result(result.vector, result.vlen); 00295 00296 if (m_df<2.0) 00297 eigen_result=CMath::INFTY*VectorXd::Ones(result.vlen); 00298 else 00299 { 00300 eigen_result+=CMath::sq(m_sigma)*m_df/(m_df-2.0)* 00301 VectorXd::Ones(result.vlen); 00302 } 00303 00304 return result; 00305 } 00306 00307 SGVector<float64_t> CStudentsTLikelihood::get_log_probability_f(const CLabels* lab, 00308 SGVector<float64_t> func) const 00309 { 00310 // check the parameters 00311 REQUIRE(lab, "Labels are required (lab should not be NULL)\n") 00312 REQUIRE(lab->get_label_type()==LT_REGRESSION, 00313 "Labels must be type of CRegressionLabels\n") 00314 REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match " 00315 "length of the function vector\n") 00316 00317 Map<VectorXd> eigen_f(func.vector, func.vlen); 00318 00319 SGVector<float64_t> r(func.vlen); 00320 Map<VectorXd> eigen_r(r.vector, r.vlen); 00321 00322 SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels(); 00323 Map<VectorXd> eigen_y(y.vector, y.vlen); 00324 00325 // compute lZ=log(gamma(df/2+1/2))-log(gamma(df/2))-log(df*pi*sigma^2)/2 00326 VectorXd eigen_lZ=(CStatistics::lgamma(m_df/2.0+0.5)- 00327 CStatistics::lgamma(m_df/2.0)-log(m_df*CMath::PI*CMath::sq(m_sigma))/2.0)* 00328 VectorXd::Ones(r.vlen); 00329 00330 // compute log probability: lp=lZ-(df+1)*log(1+(y-f).^2./(df*sigma^2))/2 00331 eigen_r=eigen_y-eigen_f; 00332 eigen_r=eigen_r.cwiseProduct(eigen_r)/(m_df*CMath::sq(m_sigma)); 00333 eigen_r=eigen_lZ-(m_df+1)* 00334 (eigen_r+VectorXd::Ones(r.vlen)).array().log().matrix()/2.0; 00335 00336 return r; 00337 } 00338 00339 SGVector<float64_t> CStudentsTLikelihood::get_log_probability_derivative_f( 00340 const CLabels* lab, SGVector<float64_t> func, index_t i) const 00341 { 00342 // check the parameters 00343 REQUIRE(lab, "Labels are required (lab should not be NULL)\n") 00344 REQUIRE(lab->get_label_type()==LT_REGRESSION, 00345 "Labels must be type of CRegressionLabels\n") 00346 REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match " 00347 "length of the function vector\n") 00348 REQUIRE(i>=1 && i<=3, "Index for derivative should be 1, 2 or 3\n") 00349 00350 Map<VectorXd> eigen_f(func.vector, func.vlen); 00351 00352 SGVector<float64_t> r(func.vlen); 00353 Map<VectorXd> eigen_r(r.vector, r.vlen); 00354 00355 SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels(); 00356 Map<VectorXd> eigen_y(y.vector, y.vlen); 00357 00358 // compute r=y-f, r2=r.^2 00359 eigen_r=eigen_y-eigen_f; 00360 VectorXd eigen_r2=eigen_r.cwiseProduct(eigen_r); 00361 00362 // compute a=(y-f).^2+df*sigma^2 00363 VectorXd a=eigen_r2+VectorXd::Ones(r.vlen)*m_df*CMath::sq(m_sigma); 00364 00365 if (i==1) 00366 { 00367 // compute first derivative of log probability wrt f: 00368 // dlp=(df+1)*(y-f)./a 00369 eigen_r=(m_df+1)*eigen_r.cwiseQuotient(a); 00370 } 00371 else if (i==2) 00372 { 00373 // compute second derivative of log probability wrt f: 00374 // d2lp=(df+1)*((y-f)^2-df*sigma^2)./a.^2 00375 VectorXd b=eigen_r2-VectorXd::Ones(r.vlen)*m_df*CMath::sq(m_sigma); 00376 00377 eigen_r=(m_df+1)*b.cwiseQuotient(a.cwiseProduct(a)); 00378 } 00379 else if (i==3) 00380 { 00381 // compute third derivative of log probability wrt f: 00382 // d3lp=(f+1)*2*(y-f).*((y-f)^2-3*f*sigma^2)./a.^3 00383 VectorXd c=eigen_r2-VectorXd::Ones(r.vlen)*3*m_df*CMath::sq(m_sigma); 00384 VectorXd a2=a.cwiseProduct(a); 00385 00386 eigen_r=(m_df+1)*2*eigen_r.cwiseProduct(c).cwiseQuotient( 00387 a2.cwiseProduct(a)); 00388 } 00389 00390 return r; 00391 } 00392 00393 SGVector<float64_t> CStudentsTLikelihood::get_first_derivative(const CLabels* lab, 00394 SGVector<float64_t> func, const TParameter* param) const 00395 { 00396 // check the parameters 00397 REQUIRE(lab, "Labels are required (lab should not be NULL)\n") 00398 REQUIRE(lab->get_label_type()==LT_REGRESSION, 00399 "Labels must be type of CRegressionLabels\n") 00400 REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match " 00401 "length of the function vector\n") 00402 00403 SGVector<float64_t> r(func.vlen); 00404 Map<VectorXd> eigen_r(r.vector, r.vlen); 00405 00406 Map<VectorXd> eigen_f(func.vector, func.vlen); 00407 00408 SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels(); 00409 Map<VectorXd> eigen_y(y.vector, y.vlen); 00410 00411 // compute r=y-f and r2=(y-f).^2 00412 eigen_r=eigen_y-eigen_f; 00413 VectorXd eigen_r2=eigen_r.cwiseProduct(eigen_r); 00414 00415 if (!strcmp(param->m_name, "df")) 00416 { 00417 // compute derivative of log probability wrt df: 00418 // lp_ddf=df*(dloggamma(df/2+1/2)-dloggamma(df/2))/2-1/2- 00419 // df*log(1+r2/(df*sigma^2))/2 +(df/2+1/2)*r2./(df*sigma^2+r2) 00420 eigen_r=(m_df*(CStatistics::dlgamma(m_df*0.5+0.5)- 00421 CStatistics::dlgamma(m_df*0.5))*0.5-0.5)*VectorXd::Ones(r.vlen); 00422 00423 eigen_r-=m_df*(VectorXd::Ones(r.vlen)+ 00424 eigen_r2/(m_df*CMath::sq(m_sigma))).array().log().matrix()/2.0; 00425 00426 eigen_r+=(m_df/2.0+0.5)*eigen_r2.cwiseQuotient( 00427 eigen_r2+VectorXd::Ones(r.vlen)*(m_df*CMath::sq(m_sigma))); 00428 00429 eigen_r*=(1.0-1.0/m_df); 00430 00431 return r; 00432 } 00433 else if (!strcmp(param->m_name, "sigma")) 00434 { 00435 // compute derivative of log probability wrt sigma: 00436 // lp_dsigma=(df+1)*r2./a-1 00437 eigen_r=(m_df+1)*eigen_r2.cwiseQuotient(eigen_r2+ 00438 VectorXd::Ones(r.vlen)*(m_df*CMath::sq(m_sigma))); 00439 eigen_r-=VectorXd::Ones(r.vlen); 00440 00441 return r; 00442 } 00443 00444 return SGVector<float64_t>(); 00445 } 00446 00447 SGVector<float64_t> CStudentsTLikelihood::get_second_derivative(const CLabels* lab, 00448 SGVector<float64_t> func, const TParameter* param) const 00449 { 00450 // check the parameters 00451 REQUIRE(lab, "Labels are required (lab should not be NULL)\n") 00452 REQUIRE(lab->get_label_type()==LT_REGRESSION, 00453 "Labels must be type of CRegressionLabels\n") 00454 REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match " 00455 "length of the function vector\n") 00456 00457 SGVector<float64_t> r(func.vlen); 00458 Map<VectorXd> eigen_r(r.vector, r.vlen); 00459 00460 Map<VectorXd> eigen_f(func.vector, func.vlen); 00461 00462 SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels(); 00463 Map<VectorXd> eigen_y(y.vector, y.vlen); 00464 00465 // compute r=y-f and r2=(y-f).^2 00466 eigen_r=eigen_y-eigen_f; 00467 VectorXd eigen_r2=eigen_r.cwiseProduct(eigen_r); 00468 00469 // compute a=r+sigma^2*df and a2=a.^2 00470 VectorXd a=eigen_r2+CMath::sq(m_sigma)*m_df*VectorXd::Ones(r.vlen); 00471 VectorXd a2=a.cwiseProduct(a); 00472 00473 if (!strcmp(param->m_name, "df")) 00474 { 00475 // compute derivative of first derivative of log probability wrt df: 00476 // dlp_ddf=df*r.*(a-sigma^2*(df+1))./a2 00477 eigen_r=m_df*eigen_r.cwiseProduct(a-CMath::sq(m_sigma)*(m_df+1.0)* 00478 VectorXd::Ones(r.vlen)).cwiseQuotient(a2); 00479 eigen_r*=(1.0-1.0/m_df); 00480 00481 return r; 00482 } 00483 else if (!strcmp(param->m_name, "sigma")) 00484 { 00485 // compute derivative of first derivative of log probability wrt sigma: 00486 // dlp_dsigma=-(df+1)*2*df*sigma^2*r./a2 00487 eigen_r=-(m_df+1.0)*2*m_df*CMath::sq(m_sigma)* 00488 eigen_r.cwiseQuotient(a2); 00489 00490 return r; 00491 } 00492 00493 return SGVector<float64_t>(); 00494 } 00495 00496 SGVector<float64_t> CStudentsTLikelihood::get_third_derivative(const CLabels* lab, 00497 SGVector<float64_t> func, const TParameter* param) const 00498 { 00499 // check the parameters 00500 REQUIRE(lab, "Labels are required (lab should not be NULL)\n") 00501 REQUIRE(lab->get_label_type()==LT_REGRESSION, 00502 "Labels must be type of CRegressionLabels\n") 00503 REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match " 00504 "length of the function vector\n") 00505 00506 SGVector<float64_t> r(func.vlen); 00507 Map<VectorXd> eigen_r(r.vector, r.vlen); 00508 00509 Map<VectorXd> eigen_f(func.vector, func.vlen); 00510 00511 SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels(); 00512 Map<VectorXd> eigen_y(y.vector, y.vlen); 00513 00514 // compute r=y-f and r2=(y-f).^2 00515 eigen_r=eigen_y-eigen_f; 00516 VectorXd eigen_r2=eigen_r.cwiseProduct(eigen_r); 00517 00518 // compute a=r+sigma^2*df and a3=a.^3 00519 VectorXd a=eigen_r2+CMath::sq(m_sigma)*m_df*VectorXd::Ones(r.vlen); 00520 VectorXd a3=(a.cwiseProduct(a)).cwiseProduct(a); 00521 00522 if (!strcmp(param->m_name, "df")) 00523 { 00524 // compute derivative of second derivative of log probability wrt df: 00525 // d2lp_ddf=df*(r2.*(r2-3*sigma^2*(1+df))+df*sigma^4)./a3 00526 float64_t sigma2=CMath::sq(m_sigma); 00527 00528 eigen_r=m_df*(eigen_r2.cwiseProduct(eigen_r2-3*sigma2*(1.0+m_df)* 00529 VectorXd::Ones(r.vlen))+(m_df*CMath::sq(sigma2))*VectorXd::Ones(r.vlen)); 00530 eigen_r=eigen_r.cwiseQuotient(a3); 00531 00532 eigen_r*=(1.0-1.0/m_df); 00533 00534 return r; 00535 } 00536 else if (!strcmp(param->m_name, "sigma")) 00537 { 00538 // compute derivative of second derivative of log probability wrt sigma: 00539 // d2lp_dsigma=(df+1)*2*df*sigma^2*(a-4*r2)./a3 00540 eigen_r=(m_df+1.0)*2*m_df*CMath::sq(m_sigma)* 00541 (a-4.0*eigen_r2).cwiseQuotient(a3); 00542 00543 return r; 00544 } 00545 00546 return SGVector<float64_t>(); 00547 } 00548 00549 SGVector<float64_t> CStudentsTLikelihood::get_log_zeroth_moments( 00550 SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const 00551 { 00552 SGVector<float64_t> y; 00553 00554 if (lab) 00555 { 00556 REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()), 00557 "Length of the vector of means (%d), length of the vector of " 00558 "variances (%d) and number of labels (%d) should be the same\n", 00559 mu.vlen, s2.vlen, lab->get_num_labels()) 00560 REQUIRE(lab->get_label_type()==LT_REGRESSION, 00561 "Labels must be type of CRegressionLabels\n") 00562 00563 y=((CRegressionLabels*)lab)->get_labels(); 00564 } 00565 else 00566 { 00567 REQUIRE(mu.vlen==s2.vlen, "Length of the vector of means (%d) and " 00568 "length of the vector of variances (%d) should be the same\n", 00569 mu.vlen, s2.vlen) 00570 00571 y=SGVector<float64_t>(mu.vlen); 00572 y.set_const(1.0); 00573 } 00574 00575 // create an object of normal pdf 00576 CNormalPDF* f=new CNormalPDF(); 00577 00578 // create an object of Student's t pdf 00579 CStudentsTPDF* g=new CStudentsTPDF(); 00580 00581 g->set_nu(m_df); 00582 g->set_sigma(m_sigma); 00583 00584 // create an object of product of Student's-t pdf and normal pdf 00585 CProductFunction* h=new CProductFunction(f, g); 00586 SG_REF(h); 00587 00588 // compute probabilities using numerical integration 00589 SGVector<float64_t> r(mu.vlen); 00590 00591 for (index_t i=0; i<mu.vlen; i++) 00592 { 00593 // set normal pdf parameters 00594 f->set_mu(mu[i]); 00595 f->set_sigma(CMath::sqrt(s2[i])); 00596 00597 // set Stundent's-t pdf parameters 00598 g->set_mu(y[i]); 00599 00600 // evaluate integral on (-inf, inf) 00601 r[i]=CIntegration::integrate_quadgk(h, -CMath::INFTY, mu[i])+ 00602 CIntegration::integrate_quadgk(h, mu[i], CMath::INFTY); 00603 } 00604 00605 SG_UNREF(h); 00606 00607 r.log(); 00608 00609 return r; 00610 } 00611 00612 float64_t CStudentsTLikelihood::get_first_moment(SGVector<float64_t> mu, 00613 SGVector<float64_t> s2, const CLabels *lab, index_t i) const 00614 { 00615 // check the parameters 00616 REQUIRE(lab, "Labels are required (lab should not be NULL)\n") 00617 REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()), 00618 "Length of the vector of means (%d), length of the vector of " 00619 "variances (%d) and number of labels (%d) should be the same\n", 00620 mu.vlen, s2.vlen, lab->get_num_labels()) 00621 REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i) 00622 REQUIRE(lab->get_label_type()==LT_REGRESSION, 00623 "Labels must be type of CRegressionLabels\n") 00624 00625 SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels(); 00626 00627 // create an object of normal pdf 00628 CNormalPDF* f=new CNormalPDF(mu[i], CMath::sqrt(s2[i])); 00629 00630 // create an object of Student's t pdf 00631 CStudentsTPDF* g=new CStudentsTPDF(m_sigma, m_df, y[i]); 00632 00633 // create an object of h(x)=N(x|mu,sigma)*t(x|mu,sigma,nu) 00634 CProductFunction* h=new CProductFunction(f, g); 00635 00636 // create an object of k(x)=x*N(x|mu,sigma)*t(x|mu,sigma,nu) 00637 CProductFunction* k=new CProductFunction(new CLinearFunction(), h); 00638 SG_REF(k); 00639 00640 // compute Z = \int N(x|mu,sigma)*t(x|mu,sigma,nu) dx 00641 float64_t Z=CIntegration::integrate_quadgk(h, -CMath::INFTY, mu[i])+ 00642 CIntegration::integrate_quadgk(h, mu[i], CMath::INFTY); 00643 00644 // compute 1st moment: 00645 // E[x] = Z^-1 * \int x*N(x|mu,sigma)*t(x|mu,sigma,nu)dx 00646 float64_t Ex=(CIntegration::integrate_quadgk(k, -CMath::INFTY, mu[i])+ 00647 CIntegration::integrate_quadgk(k, mu[i], CMath::INFTY))/Z; 00648 00649 SG_UNREF(k); 00650 00651 return Ex; 00652 } 00653 00654 float64_t CStudentsTLikelihood::get_second_moment(SGVector<float64_t> mu, 00655 SGVector<float64_t> s2, const CLabels *lab, index_t i) const 00656 { 00657 // check the parameters 00658 REQUIRE(lab, "Labels are required (lab should not be NULL)\n") 00659 REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()), 00660 "Length of the vector of means (%d), length of the vector of " 00661 "variances (%d) and number of labels (%d) should be the same\n", 00662 mu.vlen, s2.vlen, lab->get_num_labels()) 00663 REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i) 00664 REQUIRE(lab->get_label_type()==LT_REGRESSION, 00665 "Labels must be type of CRegressionLabels\n") 00666 00667 SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels(); 00668 00669 // create an object of normal pdf 00670 CNormalPDF* f=new CNormalPDF(mu[i], CMath::sqrt(s2[i])); 00671 00672 // create an object of Student's t pdf 00673 CStudentsTPDF* g=new CStudentsTPDF(m_sigma, m_df, y[i]); 00674 00675 // create an object of h(x)=N(x|mu,sigma)*t(x|mu,sigma,nu) 00676 CProductFunction* h=new CProductFunction(f, g); 00677 00678 // create an object of k(x)=x*N(x|mu,sigma)*t(x|mu,sigma,nu) 00679 CProductFunction* k=new CProductFunction(new CLinearFunction(), h); 00680 SG_REF(k); 00681 00682 // create an object of p(x)=x^2*N(x|mu,sigma^2)*t(x|mu,sigma,nu) 00683 CProductFunction* p=new CProductFunction(new CQuadraticFunction(), h); 00684 SG_REF(p); 00685 00686 // compute Z = \int N(x|mu,sigma)*t(x|mu,sigma,nu) dx 00687 float64_t Z=CIntegration::integrate_quadgk(h, -CMath::INFTY, mu[i])+ 00688 CIntegration::integrate_quadgk(h, mu[i], CMath::INFTY); 00689 00690 // compute 1st moment: 00691 // E[x] = Z^-1 * \int x*N(x|mu,sigma)*t(x|mu,sigma,nu)dx 00692 float64_t Ex=(CIntegration::integrate_quadgk(k, -CMath::INFTY, mu[i])+ 00693 CIntegration::integrate_quadgk(k, mu[i], CMath::INFTY))/Z; 00694 00695 // compute E[x^2] = Z^-1 * \int x^2*N(x|mu,sigma)*t(x|mu,sigma,nu)dx 00696 float64_t Ex2=(CIntegration::integrate_quadgk(p, -CMath::INFTY, mu[i])+ 00697 CIntegration::integrate_quadgk(p, mu[i], CMath::INFTY))/Z; 00698 00699 SG_UNREF(k); 00700 SG_UNREF(p); 00701 00702 // return 2nd moment: Var[x]=E[x^2]-E[x]^2 00703 return Ex2-CMath::sq(Ex); 00704 } 00705 } 00706 00707 #endif /* HAVE_EIGEN3 */