SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
StudentsTLikelihood.cpp
Go to the documentation of this file.
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 */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation