SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
GaussianLikelihood.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 
00012 #include <shogun/machine/gp/GaussianLikelihood.h>
00013 
00014 #ifdef HAVE_EIGEN3
00015 
00016 #include <shogun/labels/RegressionLabels.h>
00017 #include <shogun/mathematics/Math.h>
00018 #include <shogun/mathematics/eigen3.h>
00019 
00020 using namespace shogun;
00021 using namespace Eigen;
00022 
00023 CGaussianLikelihood::CGaussianLikelihood() : CLikelihoodModel()
00024 {
00025     init();
00026 }
00027 
00028 CGaussianLikelihood::CGaussianLikelihood(float64_t sigma) : CLikelihoodModel()
00029 {
00030     REQUIRE(sigma>0.0, "Standard deviation must be greater than zero\n")
00031     init();
00032     m_sigma=sigma;
00033 }
00034 
00035 void CGaussianLikelihood::init()
00036 {
00037     m_sigma=1.0;
00038     SG_ADD(&m_sigma, "sigma", "Observation noise", MS_AVAILABLE, GRADIENT_AVAILABLE);
00039 }
00040 
00041 CGaussianLikelihood::~CGaussianLikelihood()
00042 {
00043 }
00044 
00045 CGaussianLikelihood* CGaussianLikelihood::obtain_from_generic(
00046         CLikelihoodModel* lik)
00047 {
00048     ASSERT(lik!=NULL);
00049 
00050     if (lik->get_model_type()!=LT_GAUSSIAN)
00051         SG_SERROR("Provided likelihood is not of type CGaussianLikelihood!\n")
00052 
00053     SG_REF(lik);
00054     return (CGaussianLikelihood*)lik;
00055 }
00056 
00057 SGVector<float64_t> CGaussianLikelihood::get_predictive_means(
00058         SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
00059 {
00060     return SGVector<float64_t>(mu);
00061 }
00062 
00063 SGVector<float64_t> CGaussianLikelihood::get_predictive_variances(
00064         SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
00065 {
00066     SGVector<float64_t> result(s2);
00067     Map<VectorXd> eigen_result(result.vector, result.vlen);
00068 
00069     eigen_result=eigen_result.array()+CMath::sq(m_sigma);
00070 
00071     return result;
00072 }
00073 
00074 SGVector<float64_t> CGaussianLikelihood::get_log_probability_f(const CLabels* lab,
00075         SGVector<float64_t> func) const
00076 {
00077     // check the parameters
00078     REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
00079     REQUIRE(lab->get_label_type()==LT_REGRESSION,
00080             "Labels must be type of CRegressionLabels\n")
00081     REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
00082             "length of the function vector\n")
00083 
00084     Map<VectorXd> eigen_f(func.vector, func.vlen);
00085 
00086     SGVector<float64_t> result(func.vlen);
00087     Map<VectorXd> eigen_result(result.vector, result.vlen);
00088 
00089     SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels();
00090     Map<VectorXd> eigen_y(y.vector, y.vlen);
00091 
00092     // compute log probability: lp=-(y-f).^2./sigma^2/2-log(2*pi*sigma^2)/2
00093     eigen_result=eigen_y-eigen_f;
00094     eigen_result=-eigen_result.cwiseProduct(eigen_result)/(2*CMath::sq(m_sigma))-
00095         VectorXd::Ones(result.vlen)*log(2*CMath::PI*CMath::sq(m_sigma))/2.0;
00096 
00097     return result;
00098 }
00099 
00100 SGVector<float64_t> CGaussianLikelihood::get_log_probability_derivative_f(
00101         const CLabels* lab, SGVector<float64_t> func, index_t i) const
00102 {
00103     // check the parameters
00104     REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
00105     REQUIRE(lab->get_label_type()==LT_REGRESSION,
00106             "Labels must be type of CRegressionLabels\n")
00107     REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
00108             "length of the function vector\n")
00109     REQUIRE(i>=1 && i<=3, "Index for derivative should be 1, 2 or 3\n")
00110 
00111     Map<VectorXd> eigen_f(func.vector, func.vlen);
00112 
00113     SGVector<float64_t> result(func.vlen);
00114     Map<VectorXd> eigen_result(result.vector, result.vlen);
00115 
00116     SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels();
00117     Map<VectorXd> eigen_y(y.vector, y.vlen);
00118 
00119     // set result=y-f
00120     eigen_result=eigen_y-eigen_f;
00121 
00122     // compute derivatives of log probability wrt f
00123     if (i == 1)
00124         eigen_result/=CMath::sq(m_sigma);
00125     else if (i == 2)
00126         eigen_result=-VectorXd::Ones(result.vlen)/CMath::sq(m_sigma);
00127     else if (i == 3)
00128         eigen_result=VectorXd::Zero(result.vlen);
00129 
00130     return result;
00131 }
00132 
00133 SGVector<float64_t> CGaussianLikelihood::get_first_derivative(const CLabels* lab,
00134         SGVector<float64_t> func, const TParameter* param) const
00135 {
00136     REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
00137     REQUIRE(lab->get_label_type()==LT_REGRESSION,
00138             "Labels must be type of CRegressionLabels\n")
00139     REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
00140             "length of the function vector\n")
00141 
00142     Map<VectorXd> eigen_f(func.vector, func.vlen);
00143 
00144     SGVector<float64_t> result(func.vlen);
00145     Map<VectorXd> eigen_result(result.vector, result.vlen);
00146 
00147     if (strcmp(param->m_name, "sigma"))
00148         return SGVector<float64_t>();
00149 
00150     SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels();
00151     Map<VectorXd> eigen_y(y.vector, y.vlen);
00152 
00153     // compute derivative of log probability wrt sigma:
00154     // lp_dsigma=(y-f).^2/sigma^2-1
00155     eigen_result=eigen_y-eigen_f;
00156     eigen_result=eigen_result.cwiseProduct(eigen_result)/CMath::sq(m_sigma);
00157     eigen_result-=VectorXd::Ones(result.vlen);
00158 
00159     return result;
00160 }
00161 
00162 SGVector<float64_t> CGaussianLikelihood::get_second_derivative(const CLabels* lab,
00163         SGVector<float64_t> func, const TParameter* param) const
00164 {
00165     REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
00166     REQUIRE(lab->get_label_type()==LT_REGRESSION,
00167             "Labels must be type of CRegressionLabels\n")
00168     REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
00169             "length of the function vector\n")
00170 
00171     if (strcmp(param->m_name, "sigma"))
00172         return SGVector<float64_t>();
00173 
00174     Map<VectorXd> eigen_f(func.vector, func.vlen);
00175 
00176     SGVector<float64_t> result(func.vlen);
00177     Map<VectorXd> eigen_result(result.vector, result.vlen);
00178 
00179     SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels();
00180     Map<VectorXd> eigen_y(y.vector, y.vlen);
00181 
00182     // compute derivative of the first derivative of log probability wrt sigma:
00183     // dlp_dsigma=2*(f-y)/sigma^2
00184     eigen_result=2*(eigen_f-eigen_y)/CMath::sq(m_sigma);
00185 
00186     return result;
00187 }
00188 
00189 SGVector<float64_t> CGaussianLikelihood::get_third_derivative(const CLabels* lab,
00190         SGVector<float64_t> func, const TParameter* param) const
00191 {
00192     REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
00193     REQUIRE(lab->get_label_type()==LT_REGRESSION,
00194             "Labels must be type of CRegressionLabels\n")
00195     REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
00196             "length of the function vector\n")
00197 
00198     if (strcmp(param->m_name, "sigma"))
00199         return SGVector<float64_t>();
00200 
00201     Map<VectorXd> eigen_f(func.vector, func.vlen);
00202 
00203     SGVector<float64_t> result(func.vlen);
00204     Map<VectorXd> eigen_result(result.vector, result.vlen);
00205 
00206     // compute derivative of the second derivative of log probability wrt sigma:
00207     // d2lp_dsigma=1/sigma^2
00208     eigen_result=2*VectorXd::Ones(result.vlen)/CMath::sq(m_sigma);
00209 
00210     return result;
00211 }
00212 
00213 SGVector<float64_t> CGaussianLikelihood::get_log_zeroth_moments(
00214         SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels *lab) const
00215 {
00216     SGVector<float64_t> y;
00217 
00218     if (lab)
00219     {
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(lab->get_label_type()==LT_REGRESSION,
00225                 "Labels must be type of CRegressionLabels\n")
00226 
00227         y=((CRegressionLabels*)lab)->get_labels();
00228     }
00229     else
00230     {
00231         REQUIRE(mu.vlen==s2.vlen, "Length of the vector of means (%d) and "
00232                 "length of the vector of variances (%d) should be the same\n",
00233                 mu.vlen, s2.vlen)
00234 
00235         y=SGVector<float64_t>(mu.vlen);
00236         y.set_const(1.0);
00237     }
00238 
00239     // create eigen representation of y, mu and s2
00240     Map<VectorXd> eigen_mu(mu.vector, mu.vlen);
00241     Map<VectorXd> eigen_s2(s2.vector, s2.vlen);
00242     Map<VectorXd> eigen_y(y.vector, y.vlen);
00243 
00244     SGVector<float64_t> result(mu.vlen);
00245     Map<VectorXd> eigen_result(result.vector, result.vlen);
00246 
00247     // compule lZ=-(y-mu).^2./(sn2+s2)/2-log(2*pi*(sn2+s2))/2
00248     eigen_s2=eigen_s2.array()+CMath::sq(m_sigma);
00249     eigen_result=-(eigen_y-eigen_mu).array().square()/(2.0*eigen_s2.array())-
00250         (2.0*CMath::PI*eigen_s2.array()).log()/2.0;
00251 
00252     return result;
00253 }
00254 
00255 float64_t CGaussianLikelihood::get_first_moment(SGVector<float64_t> mu,
00256         SGVector<float64_t> s2, const CLabels *lab, index_t i) const
00257 {
00258     // check the parameters
00259     REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
00260     REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
00261             "Length of the vector of means (%d), length of the vector of "
00262             "variances (%d) and number of labels (%d) should be the same\n",
00263             mu.vlen, s2.vlen, lab->get_num_labels())
00264     REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
00265     REQUIRE(lab->get_label_type()==LT_REGRESSION,
00266             "Labels must be type of CRegressionLabels\n")
00267 
00268     SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels();
00269 
00270     // compute 1st moment
00271     float64_t Ex=mu[i]+s2[i]*(y[i]-mu[i])/(CMath::sq(m_sigma)+s2[i]);
00272 
00273     return Ex;
00274 }
00275 
00276 float64_t CGaussianLikelihood::get_second_moment(SGVector<float64_t> mu,
00277         SGVector<float64_t> s2, const CLabels *lab, index_t i) const
00278 {
00279     // check the parameters
00280     REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
00281     REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
00282             "Length of the vector of means (%d), length of the vector of "
00283             "variances (%d) and number of labels (%d) should be the same\n",
00284             mu.vlen, s2.vlen, lab->get_num_labels())
00285     REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
00286     REQUIRE(lab->get_label_type()==LT_REGRESSION,
00287             "Labels must be type of CRegressionLabels\n")
00288 
00289     // compute 2nd moment
00290     float64_t Var=s2[i]-CMath::sq(s2[i])/(CMath::sq(m_sigma)+s2[i]);
00291 
00292     return Var;
00293 }
00294 
00295 #endif /* HAVE_EIGEN3 */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation