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 */ 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 */