SHOGUN  6.1.3
LogitLikelihood.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (c) The Shogun Machine Learning Toolbox
3  * Written (W) 2013 Roman Votyakov
4  * All rights reserved.
5  *
6  * Redistribution and use in source and binary forms, with or without
7  * modification, are permitted provided that the following conditions are met:
8  *
9  * 1. Redistributions of source code must retain the above copyright notice, this
10  * list of conditions and the following disclaimer.
11  * 2. Redistributions in binary form must reproduce the above copyright notice,
12  * this list of conditions and the following disclaimer in the documentation
13  * and/or other materials provided with the distribution.
14  *
15  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
19  * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25  *
26  * The views and conclusions contained in the software and documentation are those
27  * of the authors and should not be interpreted as representing official policies,
28  * either expressed or implied, of the Shogun Development Team.
29  *
30  */
32 
33 
35 #ifdef USE_GPL_SHOGUN
36 #include <shogun/mathematics/Integration.h>
37 #endif //USE_GPL_SHOGUN
40 
41 using namespace shogun;
42 using namespace Eigen;
43 
44 namespace shogun
45 {
46 
47 #ifndef DOXYGEN_SHOULD_SKIP_THIS
48 
50 class CNormalPDF : public CFunction
51 {
52 public:
54  CNormalPDF()
55  {
56  m_mu=0.0;
57  m_sigma=1.0;
58  }
59 
65  CNormalPDF(float64_t mu, float64_t sigma)
66  {
67  m_mu=mu;
68  m_sigma=sigma;
69  }
70 
75  void set_mu(float64_t mu) { m_mu=mu; }
76 
81  void set_sigma(float64_t sigma) { m_sigma=sigma; }
82 
89  virtual float64_t operator() (float64_t x)
90  {
91  return (1.0/(CMath::sqrt(2*CMath::PI)*m_sigma))*
92  CMath::exp(-CMath::sq(x-m_mu)/(2.0*CMath::sq(m_sigma)));
93  }
94 
95 private:
96  /* mean */
97  float64_t m_mu;
98 
99  /* standard deviation */
100  float64_t m_sigma;
101 };
102 
104 class CSigmoidFunction : public CFunction
105 {
106 public:
108  CSigmoidFunction()
109  {
110  m_a=0.0;
111  }
112 
117  CSigmoidFunction(float64_t a)
118  {
119  m_a=a;
120  }
121 
126  void set_a(float64_t a) { m_a=a; }
127 
134  virtual float64_t operator() (float64_t x)
135  {
136  return 1.0/(1.0+CMath::exp(-m_a*x));
137  }
138 
139 private:
141  float64_t m_a;
142 };
143 
145 class CProductFunction : public CFunction
146 {
147 public:
153  CProductFunction(CFunction* f, CFunction* g)
154  {
155  SG_REF(f);
156  SG_REF(g);
157  m_f=f;
158  m_g=g;
159  }
160 
161  virtual ~CProductFunction()
162  {
163  SG_UNREF(m_f);
164  SG_UNREF(m_g);
165  }
166 
173  virtual float64_t operator() (float64_t x)
174  {
175  return (*m_f)(x)*(*m_g)(x);
176  }
177 
178 private:
180  CFunction* m_f;
182  CFunction* m_g;
183 };
184 
186 class CLinearFunction : public CFunction
187 {
188 public:
190  CLinearFunction() { }
191 
192  virtual ~CLinearFunction() { }
193 
200  virtual float64_t operator() (float64_t x)
201  {
202  return x;
203  }
204 };
205 
207 class CQuadraticFunction : public CFunction
208 {
209 public:
211  CQuadraticFunction() { }
212 
213  virtual ~CQuadraticFunction() { }
214 
221  virtual float64_t operator() (float64_t x)
222  {
223  return CMath::sq(x);
224  }
225 };
226 
227 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
228 
230 {
231 }
232 
234 {
235 }
236 
238  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
239 {
241  Map<VectorXd> eigen_lp(lp.vector, lp.vlen);
242 
244  Map<VectorXd> eigen_r(r.vector, r.vlen);
245 
246  // evaluate predictive mean: ymu=2*p-1
247  // Note that the distribution is Bernoulli distribution with p(x=1)=p and
248  // p(x=-1)=(1-p)
249  // the mean of the Bernoulli distribution is 2*p-1
250  eigen_r=2.0*eigen_lp.array().exp()-1.0;
251 
252  return r;
253 }
254 
256  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
257 {
259  Map<VectorXd> eigen_lp(lp.vector, lp.vlen);
260 
262  Map<VectorXd> eigen_r(r.vector, r.vlen);
263 
264  // evaluate predictive variance: ys2=1-(2*p-1).^2
265  // Note that the distribution is Bernoulli distribution with p(x=1)=p and
266  // p(x=-1)=(1-p)
267  // the variance of the Bernoulli distribution is 1-(2*p-1).^2
268  eigen_r=1-(2.0*eigen_lp.array().exp()-1.0).square();
269 
270  return r;
271 }
272 
274  SGVector<float64_t> func) const
275 {
276  // check the parameters
277  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
279  "Labels must be type of CBinaryLabels\n")
280  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
281  "length of the function vector\n")
282 
283  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
284  Map<VectorXd> eigen_y(y.vector, y.vlen);
285 
286  Map<VectorXd> eigen_f(func.vector, func.vlen);
287 
288  SGVector<float64_t> r(func.vlen);
289  Map<VectorXd> eigen_r(r.vector, r.vlen);
290 
291  // compute log probability: -log(1+exp(-f.*y))
292  eigen_r=-(1.0+(-eigen_y.array()*eigen_f.array()).exp()).log();
293 
294  return r;
295 }
296 
298  const CLabels* lab, SGVector<float64_t> func, index_t i) const
299 {
300  // check the parameters
301  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
303  "Labels must be type of CBinaryLabels\n")
304  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
305  "length of the function vector\n")
306  REQUIRE(i>=1 && i<=3, "Index for derivative should be 1, 2 or 3\n")
307 
308  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
309  Map<VectorXd> eigen_y(y.vector, y.vlen);
310 
311  Map<VectorXd> eigen_f(func.vector, func.vlen);
312 
313  SGVector<float64_t> r(func.vlen);
314  Map<VectorXd> eigen_r(r.vector, r.vlen);
315 
316  // compute s(f)=1./(1+exp(-f))
317  VectorXd eigen_s=(VectorXd::Ones(func.vlen)).cwiseQuotient((1.0+
318  (-eigen_f).array().exp()).matrix());
319 
320  // compute derivatives of log probability wrt f
321  if (i == 1)
322  {
323  // compute the first derivative: dlp=(y+1)/2-s(f)
324  eigen_r=(eigen_y.array()+1.0)/2.0-eigen_s.array();
325  }
326  else if (i == 2)
327  {
328  // compute the second derivative: d2lp=-s(f).*(1-s(f))
329  eigen_r=-eigen_s.array()*(1.0-eigen_s.array());
330  }
331  else if (i == 3)
332  {
333  // compute the third derivative: d2lp=-s(f).*(1-s(f)).*(1-2*s(f))
334  eigen_r=-eigen_s.array()*(1.0-eigen_s.array())*(1.0-2.0*eigen_s.array());
335  }
336  else
337  {
338  SG_ERROR("Invalid index for derivative\n")
339  }
340 
341  return r;
342 }
343 
345  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
346 {
348 
349  if (lab)
350  {
351  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
352  "Length of the vector of means (%d), length of the vector of "
353  "variances (%d) and number of labels (%d) should be the same\n",
354  mu.vlen, s2.vlen, lab->get_num_labels())
356  "Labels must be type of CBinaryLabels\n")
357 
358  y=((CBinaryLabels*)lab)->get_labels();
359  }
360  else
361  {
362  REQUIRE(mu.vlen==s2.vlen, "Length of the vector of means (%d) and "
363  "length of the vector of variances (%d) should be the same\n",
364  mu.vlen, s2.vlen)
365 
367  y.set_const(1.0);
368  }
369 
370  // create an object of normal pdf function
371  CNormalPDF* f=new CNormalPDF();
372 
373  // create an object of sigmoid function
374  CSigmoidFunction* g=new CSigmoidFunction();
375 
376  // create an object of product of sigmoid and normal pdf functions
377  CProductFunction* h=new CProductFunction(f, g);
378  SG_REF(h);
379 
380  // compute probabilities using numerical integration
382 
383  for (index_t i=0; i<mu.vlen; i++)
384  {
385  // set normal pdf parameters
386  f->set_mu(mu[i]);
387  f->set_sigma(CMath::sqrt(s2[i]));
388 
389  // set sigmoid parameters
390  g->set_a(y[i]);
391 
392  // evaluate integral on (-inf, inf)
393 #ifdef USE_GPL_SHOGUN
394  r[i]=CIntegration::integrate_quadgk(h, -CMath::INFTY, mu[i])+
395  CIntegration::integrate_quadgk(h, mu[i], CMath::INFTY);
396 #else
398 #endif //USE_GPL_SHOGUN
399  }
400 
401  SG_UNREF(h);
402 
403  for (index_t i=0; i<r.vlen; i++)
404  r[i]=CMath::log(r[i]);
405 
406  return r;
407 }
408 
410  SGVector<float64_t> s2, const CLabels *lab, index_t i) const
411 {
412  // check the parameters
413  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
414  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
415  "Length of the vector of means (%d), length of the vector of "
416  "variances (%d) and number of labels (%d) should be the same\n",
417  mu.vlen, s2.vlen, lab->get_num_labels())
418  REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
420  "Labels must be type of CBinaryLabels\n")
421 
422  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
423 
424  // create an object of f(x)=N(x|mu,sigma^2)
425  CNormalPDF* f=new CNormalPDF(mu[i], CMath::sqrt(s2[i]));
426 
427  // create an object of g(x)=sigmoid(x)
428  CSigmoidFunction* g=new CSigmoidFunction(y[i]);
429 
430  // create an object of h(x)=N(x|mu,sigma^2)*sigmoid(x)
431  CProductFunction* h=new CProductFunction(f, g);
432 
433  // create an object of l(x)=x
434  CLinearFunction* l=new CLinearFunction();
435 
436  // create an object of k(x)=x*N(x|mu,sigma^2)*sigmoid(x)
437  CProductFunction* k=new CProductFunction(l, h);
438  SG_REF(k);
439  float64_t Ex=0;
440 #ifdef USE_GPL_SHOGUN
441  // compute Z = \int N(x|mu,sigma)*sigmoid(a*x) dx
442  float64_t Z=CIntegration::integrate_quadgk(h, -CMath::INFTY, mu[i])+
443  CIntegration::integrate_quadgk(h, mu[i], CMath::INFTY);
444 
445  // compute 1st moment: E[x] = Z^-1 * \int x*N(x|mu,sigma)*sigmoid(a*x)dx
446  Ex=(CIntegration::integrate_quadgk(k, -CMath::INFTY, mu[i])+
447  CIntegration::integrate_quadgk(k, mu[i], CMath::INFTY))/Z;
448 #else
450 #endif //USE_GPL_SHOGUN
451  SG_UNREF(k);
452 
453  return Ex;
454 }
455 
457  SGVector<float64_t> s2, const CLabels *lab, index_t i) const
458 {
459  // check the parameters
460  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
461  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
462  "Length of the vector of means (%d), length of the vector of "
463  "variances (%d) and number of labels (%d) should be the same\n",
464  mu.vlen, s2.vlen, lab->get_num_labels())
465  REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
467  "Labels must be type of CBinaryLabels\n")
468 
469  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
470 
471  // create an object of f(x)=N(x|mu,sigma^2)
472  CNormalPDF* f=new CNormalPDF(mu[i], CMath::sqrt(s2[i]));
473 
474  // create an object of g(x)=sigmoid(a*x)
475  CSigmoidFunction* g=new CSigmoidFunction(y[i]);
476 
477  // create an object of h(x)=N(x|mu,sigma^2)*sigmoid(a*x)
478  CProductFunction* h=new CProductFunction(f, g);
479 
480  // create an object of l(x)=x
481  CLinearFunction* l=new CLinearFunction();
482 
483  // create an object of k(x)=x*N(x|mu,sigma^2)*sigmoid(a*x)
484  CProductFunction* k=new CProductFunction(l, h);
485  SG_REF(k);
486 
487  // create an object of q(x)=x^2
488  CQuadraticFunction* q=new CQuadraticFunction();
489 
490  // create an object of p(x)=x^2*N(x|mu,sigma^2)*sigmoid(x)
491  CProductFunction* p=new CProductFunction(q, h);
492  SG_REF(p);
493 
494  float64_t Ex=0;
495  float64_t Ex2=0;
496 #ifdef USE_GPL_SHOGUN
497  // compute Z = \int N(x|mu,sigma)*sigmoid(a*x) dx
498  float64_t Z=CIntegration::integrate_quadgk(h, -CMath::INFTY, mu[i])+
499  CIntegration::integrate_quadgk(h, mu[i], CMath::INFTY);
500 
501  // compute 1st moment: E[x] = Z^-1 * \int x*N(x|mu,sigma)*sigmoid(a*x)dx
502  Ex=(CIntegration::integrate_quadgk(k, -CMath::INFTY, mu[i])+
503  CIntegration::integrate_quadgk(k, mu[i], CMath::INFTY))/Z;
504 
505  // compute E[x^2] = Z^-1 * \int x^2*N(x|mu,sigma)*sigmoid(a*x)dx
506  Ex2=(CIntegration::integrate_quadgk(p, -CMath::INFTY, mu[i])+
507  CIntegration::integrate_quadgk(p, mu[i], CMath::INFTY))/Z;
508 #else
510 #endif //USE_GPL_SHOGUN
511  SG_UNREF(k);
512  SG_UNREF(p);
513 
514  // return 2nd moment: Var[x]=E[x^2]-E[x]^2
515  return Ex2-CMath::sq(Ex);;
516 }
517 }
518 
virtual SGVector< float64_t > get_predictive_means(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab=NULL) const
virtual SGVector< float64_t > get_predictive_variances(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab=NULL) const
virtual ELabelType get_label_type() const =0
binary labels +1/-1
Definition: LabelTypes.h:18
int32_t index_t
Definition: common.h:72
The class Labels models labels, i.e. class assignments of objects.
Definition: Labels.h:43
static const float64_t INFTY
infinity
Definition: Math.h:1868
virtual int32_t get_num_labels() const =0
static T sqrt(T x)
Definition: Math.h:428
static T sq(T x)
Definition: Math.h:418
Definition: SGMatrix.h:25
#define SG_ERROR(...)
Definition: SGIO.h:128
#define REQUIRE(x,...)
Definition: SGIO.h:181
#define SG_REF(x)
Definition: SGObject.h:52
virtual float64_t get_second_moment(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab, index_t i) const
Class of a function of one variable.
Definition: Function.h:22
virtual SGVector< float64_t > get_log_probability_derivative_f(const CLabels *lab, SGVector< float64_t > func, index_t i) const
virtual SGVector< float64_t > get_log_zeroth_moments(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab) const
#define SG_GPL_ONLY
Definition: SGIO.h:139
double float64_t
Definition: common.h:60
void set_const(T const_elem)
Definition: SGVector.cpp:199
#define SG_UNREF(x)
Definition: SGObject.h:53
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
static float64_t exp(float64_t x)
Definition: Math.h:551
static float64_t log(float64_t v)
Definition: Math.h:714
virtual float64_t get_first_moment(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab, index_t i) const
Binary Labels for binary classification.
Definition: BinaryLabels.h:37
virtual SGVector< float64_t > get_log_probability_f(const CLabels *lab, SGVector< float64_t > func) const
The Likelihood model base class.
index_t vlen
Definition: SGVector.h:571
static const float64_t PI
Definition: Math.h:1875

SHOGUN Machine Learning Toolbox - Documentation