SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
LogitLikelihood.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  */
00009 
00010 #include <shogun/machine/gp/LogitLikelihood.h>
00011 
00012 #ifdef HAVE_EIGEN3
00013 
00014 #include <shogun/mathematics/Function.h>
00015 #include <shogun/mathematics/Integration.h>
00016 #include <shogun/labels/BinaryLabels.h>
00017 #include <shogun/mathematics/eigen3.h>
00018 
00019 using namespace shogun;
00020 using namespace Eigen;
00021 
00022 namespace shogun
00023 {
00024 
00025 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00026 
00028 class CNormalPDF : public CFunction
00029 {
00030 public:
00032     CNormalPDF()
00033     {
00034         m_mu=0.0;
00035         m_sigma=1.0;
00036     }
00037 
00043     CNormalPDF(float64_t mu, float64_t sigma)
00044     {
00045         m_mu=mu;
00046         m_sigma=sigma;
00047     }
00048 
00053     void set_mu(float64_t mu) { m_mu=mu; }
00054 
00059     void set_sigma(float64_t sigma) { m_sigma=sigma; }
00060 
00067     virtual float64_t operator() (float64_t x)
00068     {
00069         return (1.0/(CMath::sqrt(2*CMath::PI)*m_sigma))*
00070             CMath::exp(-CMath::sq(x-m_mu)/(2.0*CMath::sq(m_sigma)));
00071     }
00072 
00073 private:
00074     /* mean */
00075     float64_t m_mu;
00076 
00077     /* standard deviation */
00078     float64_t m_sigma;
00079 };
00080 
00082 class CSigmoidFunction : public CFunction
00083 {
00084 public:
00086     CSigmoidFunction()
00087     {
00088         m_a=0.0;
00089     }
00090 
00095     CSigmoidFunction(float64_t a)
00096     {
00097         m_a=a;
00098     }
00099 
00104     void set_a(float64_t a) { m_a=a; }
00105 
00112     virtual float64_t operator() (float64_t x)
00113     {
00114         return 1.0/(1.0+CMath::exp(-m_a*x));
00115     }
00116 
00117 private:
00119     float64_t m_a;
00120 };
00121 
00123 class CProductFunction : public CFunction
00124 {
00125 public:
00131     CProductFunction(CFunction* f, CFunction* g)
00132     {
00133         SG_REF(f);
00134         SG_REF(g);
00135         m_f=f;
00136         m_g=g;
00137     }
00138 
00139     virtual ~CProductFunction()
00140     {
00141         SG_UNREF(m_f);
00142         SG_UNREF(m_g);
00143     }
00144 
00151     virtual float64_t operator() (float64_t x)
00152     {
00153         return (*m_f)(x)*(*m_g)(x);
00154     }
00155 
00156 private:
00158     CFunction* m_f;
00160     CFunction* m_g;
00161 };
00162 
00164 class CLinearFunction : public CFunction
00165 {
00166 public:
00168     CLinearFunction() { }
00169 
00170     virtual ~CLinearFunction() { }
00171 
00178     virtual float64_t operator() (float64_t x)
00179     {
00180         return x;
00181     }
00182 };
00183 
00185 class CQuadraticFunction : public CFunction
00186 {
00187 public:
00189     CQuadraticFunction() { }
00190 
00191     virtual ~CQuadraticFunction() { }
00192 
00199     virtual float64_t operator() (float64_t x)
00200     {
00201         return CMath::sq(x);
00202     }
00203 };
00204 
00205 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00206 
00207 CLogitLikelihood::CLogitLikelihood() : CLikelihoodModel()
00208 {
00209 }
00210 
00211 CLogitLikelihood::~CLogitLikelihood()
00212 {
00213 }
00214 
00215 SGVector<float64_t> CLogitLikelihood::get_predictive_means(
00216         SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
00217 {
00218     SGVector<float64_t> lp=get_log_zeroth_moments(mu, s2, lab);
00219     Map<VectorXd> eigen_lp(lp.vector, lp.vlen);
00220 
00221     SGVector<float64_t> r(lp.vlen);
00222     Map<VectorXd> eigen_r(r.vector, r.vlen);
00223 
00224     // evaluate predictive mean: ymu=2*p-1
00225     eigen_r=2.0*eigen_lp.array().exp()-1.0;
00226 
00227     return r;
00228 }
00229 
00230 SGVector<float64_t> CLogitLikelihood::get_predictive_variances(
00231         SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
00232 {
00233     SGVector<float64_t> lp=get_log_zeroth_moments(mu, s2, lab);
00234     Map<VectorXd> eigen_lp(lp.vector, lp.vlen);
00235 
00236     SGVector<float64_t> r(lp.vlen);
00237     Map<VectorXd> eigen_r(r.vector, r.vlen);
00238 
00239     // evaluate predictive variance: ys2=1-(2*p-1).^2
00240     eigen_r=1-(2.0*eigen_lp.array().exp()-1.0).square();
00241 
00242     return r;
00243 }
00244 
00245 SGVector<float64_t> CLogitLikelihood::get_log_probability_f(const CLabels* lab,
00246         SGVector<float64_t> func) const
00247 {
00248     // check the parameters
00249     REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
00250     REQUIRE(lab->get_label_type()==LT_BINARY,
00251             "Labels must be type of CBinaryLabels\n")
00252     REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
00253             "length of the function vector\n")
00254 
00255     SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
00256     Map<VectorXd> eigen_y(y.vector, y.vlen);
00257 
00258     Map<VectorXd> eigen_f(func.vector, func.vlen);
00259 
00260     SGVector<float64_t> r(func.vlen);
00261     Map<VectorXd> eigen_r(r.vector, r.vlen);
00262 
00263     // compute log probability: -log(1+exp(-f.*y))
00264     eigen_r=-(1.0+(-eigen_y.array()*eigen_f.array()).exp()).log();
00265 
00266     return r;
00267 }
00268 
00269 SGVector<float64_t> CLogitLikelihood::get_log_probability_derivative_f(
00270         const CLabels* lab, SGVector<float64_t> func, index_t i) const
00271 {
00272     // check the parameters
00273     REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
00274     REQUIRE(lab->get_label_type()==LT_BINARY,
00275             "Labels must be type of CBinaryLabels\n")
00276     REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
00277             "length of the function vector\n")
00278     REQUIRE(i>=1 && i<=3, "Index for derivative should be 1, 2 or 3\n")
00279 
00280     SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
00281     Map<VectorXd> eigen_y(y.vector, y.vlen);
00282 
00283     Map<VectorXd> eigen_f(func.vector, func.vlen);
00284 
00285     SGVector<float64_t> r(func.vlen);
00286     Map<VectorXd> eigen_r(r.vector, r.vlen);
00287 
00288     // compute s(f)=1./(1+exp(-f))
00289     VectorXd eigen_s=(VectorXd::Ones(func.vlen)).cwiseQuotient((1.0+
00290         (-eigen_f).array().exp()).matrix());
00291 
00292     // compute derivatives of log probability wrt f
00293     if (i == 1)
00294     {
00295         // compute the first derivative: dlp=(y+1)/2-s(f)
00296         eigen_r=(eigen_y.array()+1.0)/2.0-eigen_s.array();
00297     }
00298     else if (i == 2)
00299     {
00300         // compute the second derivative: d2lp=-s(f).*(1-s(f))
00301         eigen_r=-eigen_s.array()*(1.0-eigen_s.array());
00302     }
00303     else if (i == 3)
00304     {
00305         // compute the third derivative: d2lp=-s(f).*(1-s(f)).*(1-2*s(f))
00306         eigen_r=-eigen_s.array()*(1.0-eigen_s.array())*(1.0-2*eigen_s.array());
00307     }
00308     else
00309     {
00310         SG_ERROR("Invalid index for derivative\n")
00311     }
00312 
00313     return r;
00314 }
00315 
00316 SGVector<float64_t> CLogitLikelihood::get_log_zeroth_moments(
00317         SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
00318 {
00319     SGVector<float64_t> y;
00320 
00321     if (lab)
00322     {
00323         REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
00324                 "Length of the vector of means (%d), length of the vector of "
00325                 "variances (%d) and number of labels (%d) should be the same\n",
00326                 mu.vlen, s2.vlen, lab->get_num_labels())
00327         REQUIRE(lab->get_label_type()==LT_BINARY,
00328                 "Labels must be type of CBinaryLabels\n")
00329 
00330         y=((CBinaryLabels*)lab)->get_labels();
00331     }
00332     else
00333     {
00334         REQUIRE(mu.vlen==s2.vlen, "Length of the vector of means (%d) and "
00335                 "length of the vector of variances (%d) should be the same\n",
00336                 mu.vlen, s2.vlen)
00337 
00338         y=SGVector<float64_t>(mu.vlen);
00339         y.set_const(1.0);
00340     }
00341 
00342     // create an object of normal pdf function
00343     CNormalPDF* f=new CNormalPDF();
00344 
00345     // create an object of sigmoid function
00346     CSigmoidFunction* g=new CSigmoidFunction();
00347 
00348     // create an object of product of sigmoid and normal pdf functions
00349     CProductFunction* h=new CProductFunction(f, g);
00350     SG_REF(h);
00351 
00352     // compute probabilities using numerical integration
00353     SGVector<float64_t> r(mu.vlen);
00354 
00355     for (index_t i=0; i<mu.vlen; i++)
00356     {
00357         // set normal pdf parameters
00358         f->set_mu(mu[i]);
00359         f->set_sigma(CMath::sqrt(s2[i]));
00360 
00361         // set sigmoid parameters
00362         g->set_a(y[i]);
00363 
00364         // evaluate integral on (-inf, inf)
00365         r[i]=CIntegration::integrate_quadgk(h, -CMath::INFTY, mu[i])+
00366             CIntegration::integrate_quadgk(h, mu[i], CMath::INFTY);
00367     }
00368 
00369     SG_UNREF(h);
00370 
00371     r.log();
00372 
00373     return r;
00374 }
00375 
00376 float64_t CLogitLikelihood::get_first_moment(SGVector<float64_t> mu,
00377         SGVector<float64_t> s2, const CLabels *lab, index_t i) const
00378 {
00379     // check the parameters
00380     REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
00381     REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
00382             "Length of the vector of means (%d), length of the vector of "
00383             "variances (%d) and number of labels (%d) should be the same\n",
00384             mu.vlen, s2.vlen, lab->get_num_labels())
00385     REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
00386     REQUIRE(lab->get_label_type()==LT_BINARY,
00387             "Labels must be type of CBinaryLabels\n")
00388 
00389     SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
00390 
00391     // create an object of f(x)=N(x|mu,sigma^2)
00392     CNormalPDF* f=new CNormalPDF(mu[i], CMath::sqrt(s2[i]));
00393 
00394     // create an object of g(x)=sigmoid(x)
00395     CSigmoidFunction* g=new CSigmoidFunction(y[i]);
00396 
00397     // create an object of h(x)=N(x|mu,sigma^2)*sigmoid(x)
00398     CProductFunction* h=new CProductFunction(f, g);
00399 
00400     // create an object of l(x)=x
00401     CLinearFunction* l=new CLinearFunction();
00402 
00403     // create an object of k(x)=x*N(x|mu,sigma^2)*sigmoid(x)
00404     CProductFunction* k=new CProductFunction(l, h);
00405     SG_REF(k);
00406 
00407     // compute Z = \int N(x|mu,sigma)*sigmoid(a*x) dx
00408     float64_t Z=CIntegration::integrate_quadgk(h, -CMath::INFTY, mu[i])+
00409         CIntegration::integrate_quadgk(h, mu[i], CMath::INFTY);
00410 
00411     // compute 1st moment: E[x] = Z^-1 * \int x*N(x|mu,sigma)*sigmoid(a*x)dx
00412     float64_t Ex=(CIntegration::integrate_quadgk(k, -CMath::INFTY, mu[i])+
00413             CIntegration::integrate_quadgk(k, mu[i], CMath::INFTY))/Z;
00414 
00415     SG_UNREF(k);
00416 
00417     return Ex;
00418 }
00419 
00420 float64_t CLogitLikelihood::get_second_moment(SGVector<float64_t> mu,
00421         SGVector<float64_t> s2, const CLabels *lab, index_t i) const
00422 {
00423     // check the parameters
00424     REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
00425     REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
00426             "Length of the vector of means (%d), length of the vector of "
00427             "variances (%d) and number of labels (%d) should be the same\n",
00428             mu.vlen, s2.vlen, lab->get_num_labels())
00429     REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
00430     REQUIRE(lab->get_label_type()==LT_BINARY,
00431             "Labels must be type of CBinaryLabels\n")
00432 
00433     SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
00434 
00435     // create an object of f(x)=N(x|mu,sigma^2)
00436     CNormalPDF* f=new CNormalPDF(mu[i], CMath::sqrt(s2[i]));
00437 
00438     // create an object of g(x)=sigmoid(a*x)
00439     CSigmoidFunction* g=new CSigmoidFunction(y[i]);
00440 
00441     // create an object of h(x)=N(x|mu,sigma^2)*sigmoid(a*x)
00442     CProductFunction* h=new CProductFunction(f, g);
00443 
00444     // create an object of l(x)=x
00445     CLinearFunction* l=new CLinearFunction();
00446 
00447     // create an object of k(x)=x*N(x|mu,sigma^2)*sigmoid(a*x)
00448     CProductFunction* k=new CProductFunction(l, h);
00449     SG_REF(k);
00450 
00451     // create an object of q(x)=x^2
00452     CQuadraticFunction* q=new CQuadraticFunction();
00453 
00454     // create an object of p(x)=x^2*N(x|mu,sigma^2)*sigmoid(x)
00455     CProductFunction* p=new CProductFunction(q, h);
00456     SG_REF(p);
00457 
00458     // compute Z = \int N(x|mu,sigma)*sigmoid(a*x) dx
00459     float64_t Z=CIntegration::integrate_quadgk(h, -CMath::INFTY, mu[i])+
00460         CIntegration::integrate_quadgk(h, mu[i], CMath::INFTY);
00461 
00462     // compute 1st moment: E[x] = Z^-1 * \int x*N(x|mu,sigma)*sigmoid(a*x)dx
00463     float64_t Ex=(CIntegration::integrate_quadgk(k, -CMath::INFTY, mu[i])+
00464             CIntegration::integrate_quadgk(k, mu[i], CMath::INFTY))/Z;
00465 
00466     // compute E[x^2] = Z^-1 * \int x^2*N(x|mu,sigma)*sigmoid(a*x)dx
00467     float64_t Ex2=(CIntegration::integrate_quadgk(p, -CMath::INFTY, mu[i])+
00468             CIntegration::integrate_quadgk(p, mu[i], CMath::INFTY))/Z;
00469 
00470     SG_UNREF(k);
00471     SG_UNREF(p);
00472 
00473     // return 2nd moment: Var[x]=E[x^2]-E[x]^2
00474     return Ex2-CMath::sq(Ex);;
00475 }
00476 }
00477 
00478 #endif /* HAVE_EIGEN3 */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation