SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
ProbitLikelihood.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/ProbitLikelihood.h>
00011 
00012 #ifdef HAVE_EIGEN3
00013 
00014 #include <shogun/labels/BinaryLabels.h>
00015 #include <shogun/mathematics/eigen3.h>
00016 #include <shogun/mathematics/Statistics.h>
00017 
00018 using namespace shogun;
00019 using namespace Eigen;
00020 
00021 CProbitLikelihood::CProbitLikelihood()
00022 {
00023 }
00024 
00025 CProbitLikelihood::~CProbitLikelihood()
00026 {
00027 }
00028 
00029 SGVector<float64_t> CProbitLikelihood::get_predictive_means(
00030         SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
00031 {
00032     SGVector<float64_t> lp=get_log_zeroth_moments(mu, s2, lab);
00033     Map<VectorXd> eigen_lp(lp.vector, lp.vlen);
00034 
00035     SGVector<float64_t> r(lp.vlen);
00036     Map<VectorXd> eigen_r(r.vector, r.vlen);
00037 
00038     // evaluate predictive mean: ymu=2*p-1
00039     eigen_r=2.0*eigen_lp.array().exp()-1.0;
00040 
00041     return r;
00042 }
00043 
00044 SGVector<float64_t> CProbitLikelihood::get_predictive_variances(
00045         SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
00046 {
00047     SGVector<float64_t> lp=get_log_zeroth_moments(mu, s2, lab);
00048     Map<VectorXd> eigen_lp(lp.vector, lp.vlen);
00049 
00050     SGVector<float64_t> r(lp.vlen);
00051     Map<VectorXd> eigen_r(r.vector, r.vlen);
00052 
00053     // evaluate predictive variance: ys2=1-(2*p-1).^2
00054     eigen_r=1-(2.0*eigen_lp.array().exp()-1.0).square();
00055 
00056     return r;
00057 }
00058 
00059 SGVector<float64_t> CProbitLikelihood::get_log_probability_f(const CLabels* lab,
00060         SGVector<float64_t> func) const
00061 {
00062     // check the parameters
00063     REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
00064     REQUIRE(lab->get_label_type()==LT_BINARY,
00065             "Labels must be type of CBinaryLabels\n")
00066     REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
00067             "length of the function vector\n")
00068 
00069     SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
00070     Map<VectorXd> eigen_y(y.vector, y.vlen);
00071 
00072     Map<VectorXd> eigen_f(func.vector, func.vlen);
00073 
00074     SGVector<float64_t> r(func.vlen);
00075     Map<VectorXd> eigen_r(r.vector, r.vlen);
00076 
00077     // compute log pobability: log(normal_cdf(f.*y))
00078     eigen_r=eigen_y.cwiseProduct(eigen_f);
00079 
00080     for (index_t i=0; i<eigen_r.size(); i++)
00081         eigen_r[i]=CStatistics::lnormal_cdf(eigen_r[i]);
00082 
00083     return r;
00084 }
00085 
00086 SGVector<float64_t> CProbitLikelihood::get_log_probability_derivative_f(
00087         const CLabels* lab, SGVector<float64_t> func, index_t i) const
00088 {
00089     // check the parameters
00090     REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
00091     REQUIRE(lab->get_label_type()==LT_BINARY,
00092             "Labels must be type of CBinaryLabels\n")
00093     REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
00094             "length of the function vector\n")
00095     REQUIRE(i>=1 && i<=3, "Index for derivative should be 1, 2 or 3\n")
00096 
00097     SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
00098     Map<VectorXd> eigen_y(y.vector, y.vlen);
00099 
00100     Map<VectorXd> eigen_f(func.vector, func.vlen);
00101 
00102     SGVector<float64_t> r(func.vlen);
00103     Map<VectorXd> eigen_r(r.vector, r.vlen);
00104 
00105     // compute ncdf=normal_cdf(y.*f)
00106     VectorXd eigen_ncdf=eigen_y.cwiseProduct(eigen_f);
00107 
00108     for (index_t j=0; j<eigen_ncdf.size(); j++)
00109         eigen_ncdf[j]=CStatistics::normal_cdf(eigen_ncdf[j]);
00110 
00111     // compute npdf=normal_pdf(f)=(1/sqrt(2*pi))*exp(-f.^2/2)
00112     VectorXd eigen_npdf=(1.0/CMath::sqrt(2.0*CMath::PI))*
00113         (-0.5*eigen_f.array().square()).exp();
00114 
00115     // compute z=npdf/ncdf
00116     VectorXd eigen_z=eigen_npdf.cwiseQuotient(eigen_ncdf);
00117 
00118     // compute derivatives of log probability wrt f
00119     if (i == 1)
00120     {
00121         // compute the first derivative: dlp=y*z
00122         eigen_r=eigen_y.cwiseProduct(eigen_z);
00123     }
00124     else if (i == 2)
00125     {
00126         // compute the second derivative: d2lp=-z.^2-y.*f.*z
00127         eigen_r=-eigen_z.array().square()-eigen_y.array()*eigen_f.array()*
00128             eigen_z.array();
00129     }
00130     else if (i == 3)
00131     {
00132         VectorXd eigen_z2=eigen_z.cwiseProduct(eigen_z);
00133         VectorXd eigen_z3=eigen_z2.cwiseProduct(eigen_z);
00134 
00135         // compute the third derivative: d3lp=2*y.*z.^3+3*f.*z.^2+z.*y.*(f.^2-1)
00136         eigen_r=2.0*eigen_y.array()*eigen_z3.array()+3.0*eigen_f.array()*
00137             eigen_z2.array()+eigen_z.array()*eigen_y.array()*
00138             (eigen_f.array().square()-1.0);
00139     }
00140 
00141     return r;
00142 }
00143 
00144 SGVector<float64_t> CProbitLikelihood::get_log_zeroth_moments(
00145         SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
00146 {
00147     SGVector<float64_t> y;
00148 
00149     if (lab)
00150     {
00151         REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
00152                 "Length of the vector of means (%d), length of the vector of "
00153                 "variances (%d) and number of labels (%d) should be the same\n",
00154                 mu.vlen, s2.vlen, lab->get_num_labels())
00155         REQUIRE(lab->get_label_type()==LT_BINARY,
00156                 "Labels must be type of CBinaryLabels\n")
00157 
00158         y=((CBinaryLabels*)lab)->get_labels();
00159     }
00160     else
00161     {
00162         REQUIRE(mu.vlen==s2.vlen, "Length of the vector of means (%d) and "
00163                 "length of the vector of variances (%d) should be the same\n",
00164                 mu.vlen, s2.vlen)
00165 
00166         y=SGVector<float64_t>(mu.vlen);
00167         y.set_const(1.0);
00168     }
00169 
00170     Map<VectorXd> eigen_y(y.vector, y.vlen);
00171     Map<VectorXd> eigen_mu(mu.vector, mu.vlen);
00172     Map<VectorXd> eigen_s2(s2.vector, s2.vlen);
00173 
00174     SGVector<float64_t> r(y.vlen);
00175     Map<VectorXd> eigen_r(r.vector, r.vlen);
00176 
00177     // compute: lp=log(normal_cdf((mu.*y)./sqrt(1+sigma^2)))
00178     eigen_r=eigen_mu.array()*eigen_y.array()/((1.0+eigen_s2.array()).sqrt());
00179 
00180     for (index_t i=0; i<eigen_r.size(); i++)
00181         eigen_r[i]=CStatistics::lnormal_cdf(eigen_r[i]);
00182 
00183     return r;
00184 }
00185 
00186 float64_t CProbitLikelihood::get_first_moment(SGVector<float64_t> mu,
00187         SGVector<float64_t> s2, const CLabels *lab, index_t i) const
00188 {
00189     // check the parameters
00190     REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
00191     REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
00192             "Length of the vector of means (%d), length of the vector of "
00193             "variances (%d) and number of labels (%d) should be the same\n",
00194             mu.vlen, s2.vlen, lab->get_num_labels())
00195     REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
00196     REQUIRE(lab->get_label_type()==LT_BINARY,
00197             "Labels must be type of CBinaryLabels\n")
00198 
00199     SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
00200 
00201     float64_t z=y[i]*mu[i]/CMath::sqrt(1.0+s2[i]);
00202 
00203     // compute ncdf=normal_cdf(z)
00204     float64_t ncdf=CStatistics::normal_cdf(z);
00205 
00206     // compute npdf=normal_pdf(z)=(1/sqrt(2*pi))*exp(-z.^2/2)
00207     float64_t npdf=(1.0/CMath::sqrt(2.0*CMath::PI))*CMath::exp(-0.5*CMath::sq(z));
00208 
00209     // compute the 1st moment: E[x] = mu + (y*s2*N(z))/(Phi(z)*sqrt(1+s2))
00210     float64_t Ex=mu[i]+(npdf/ncdf)*(y[i]*s2[i])/CMath::sqrt(1.0+s2[i]);
00211 
00212     return Ex;
00213 }
00214 
00215 float64_t CProbitLikelihood::get_second_moment(SGVector<float64_t> mu,
00216         SGVector<float64_t> s2, const CLabels *lab, index_t i) const
00217 {
00218     // check the parameters
00219     REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
00220     REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
00221             "Length of the vector of means (%d), length of the vector of "
00222             "variances (%d) and number of labels (%d) should be the same\n",
00223             mu.vlen, s2.vlen, lab->get_num_labels())
00224     REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
00225     REQUIRE(lab->get_label_type()==LT_BINARY,
00226             "Labels must be type of CBinaryLabels\n")
00227 
00228     SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
00229 
00230     float64_t z=y[i]*mu[i]/CMath::sqrt(1.0+s2[i]);
00231 
00232     // compute ncdf=normal_cdf(z)
00233     float64_t ncdf=CStatistics::normal_cdf(z);
00234 
00235     // compute npdf=normal_pdf(z)=(1/sqrt(2*pi))*exp(-z.^2/2)
00236     float64_t npdf=(1.0/CMath::sqrt(2.0*CMath::PI))*CMath::exp(-0.5*CMath::sq(z));
00237 
00238     SGVector<float64_t> r(y.vlen);
00239     Map<VectorXd> eigen_r(r.vector, r.vlen);
00240 
00241     // compute the 2nd moment:
00242     // Var[x] = s2 - (s2^2*N(z))/((1+s2)*Phi(z))*(z+N(z)/Phi(z))
00243     float64_t Var=s2[i]-(CMath::sq(s2[i])/(1.0+s2[i]))*(npdf/ncdf)*(z+(npdf/ncdf));
00244 
00245     return Var;
00246 }
00247 
00248 #endif /* HAVE_EIGEN3 */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation