SHOGUN  4.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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 
38 
39 using namespace shogun;
40 using namespace Eigen;
41 
42 namespace shogun
43 {
44 
45 #ifndef DOXYGEN_SHOULD_SKIP_THIS
46 
48 class CNormalPDF : public CFunction
49 {
50 public:
52  CNormalPDF()
53  {
54  m_mu=0.0;
55  m_sigma=1.0;
56  }
57 
63  CNormalPDF(float64_t mu, float64_t sigma)
64  {
65  m_mu=mu;
66  m_sigma=sigma;
67  }
68 
73  void set_mu(float64_t mu) { m_mu=mu; }
74 
79  void set_sigma(float64_t sigma) { m_sigma=sigma; }
80 
87  virtual float64_t operator() (float64_t x)
88  {
89  return (1.0/(CMath::sqrt(2*CMath::PI)*m_sigma))*
90  CMath::exp(-CMath::sq(x-m_mu)/(2.0*CMath::sq(m_sigma)));
91  }
92 
93 private:
94  /* mean */
95  float64_t m_mu;
96 
97  /* standard deviation */
98  float64_t m_sigma;
99 };
100 
102 class CSigmoidFunction : public CFunction
103 {
104 public:
106  CSigmoidFunction()
107  {
108  m_a=0.0;
109  }
110 
115  CSigmoidFunction(float64_t a)
116  {
117  m_a=a;
118  }
119 
124  void set_a(float64_t a) { m_a=a; }
125 
132  virtual float64_t operator() (float64_t x)
133  {
134  return 1.0/(1.0+CMath::exp(-m_a*x));
135  }
136 
137 private:
139  float64_t m_a;
140 };
141 
143 class CProductFunction : public CFunction
144 {
145 public:
151  CProductFunction(CFunction* f, CFunction* g)
152  {
153  SG_REF(f);
154  SG_REF(g);
155  m_f=f;
156  m_g=g;
157  }
158 
159  virtual ~CProductFunction()
160  {
161  SG_UNREF(m_f);
162  SG_UNREF(m_g);
163  }
164 
171  virtual float64_t operator() (float64_t x)
172  {
173  return (*m_f)(x)*(*m_g)(x);
174  }
175 
176 private:
178  CFunction* m_f;
180  CFunction* m_g;
181 };
182 
184 class CLinearFunction : public CFunction
185 {
186 public:
188  CLinearFunction() { }
189 
190  virtual ~CLinearFunction() { }
191 
198  virtual float64_t operator() (float64_t x)
199  {
200  return x;
201  }
202 };
203 
205 class CQuadraticFunction : public CFunction
206 {
207 public:
209  CQuadraticFunction() { }
210 
211  virtual ~CQuadraticFunction() { }
212 
219  virtual float64_t operator() (float64_t x)
220  {
221  return CMath::sq(x);
222  }
223 };
224 
225 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
226 
228 {
229 }
230 
232 {
233 }
234 
236  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
237 {
239  Map<VectorXd> eigen_lp(lp.vector, lp.vlen);
240 
242  Map<VectorXd> eigen_r(r.vector, r.vlen);
243 
244  // evaluate predictive mean: ymu=2*p-1
245  // Note that the distribution is Bernoulli distribution with p(x=1)=p and
246  // p(x=-1)=(1-p)
247  // the mean of the Bernoulli distribution is 2*p-1
248  eigen_r=2.0*eigen_lp.array().exp()-1.0;
249 
250  return r;
251 }
252 
254  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
255 {
257  Map<VectorXd> eigen_lp(lp.vector, lp.vlen);
258 
260  Map<VectorXd> eigen_r(r.vector, r.vlen);
261 
262  // evaluate predictive variance: ys2=1-(2*p-1).^2
263  // Note that the distribution is Bernoulli distribution with p(x=1)=p and
264  // p(x=-1)=(1-p)
265  // the variance of the Bernoulli distribution is 1-(2*p-1).^2
266  eigen_r=1-(2.0*eigen_lp.array().exp()-1.0).square();
267 
268  return r;
269 }
270 
272  SGVector<float64_t> func) const
273 {
274  // check the parameters
275  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
277  "Labels must be type of CBinaryLabels\n")
278  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
279  "length of the function vector\n")
280 
281  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
282  Map<VectorXd> eigen_y(y.vector, y.vlen);
283 
284  Map<VectorXd> eigen_f(func.vector, func.vlen);
285 
286  SGVector<float64_t> r(func.vlen);
287  Map<VectorXd> eigen_r(r.vector, r.vlen);
288 
289  // compute log probability: -log(1+exp(-f.*y))
290  eigen_r=-(1.0+(-eigen_y.array()*eigen_f.array()).exp()).log();
291 
292  return r;
293 }
294 
296  const CLabels* lab, SGVector<float64_t> func, index_t i) const
297 {
298  // check the parameters
299  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
301  "Labels must be type of CBinaryLabels\n")
302  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
303  "length of the function vector\n")
304  REQUIRE(i>=1 && i<=3, "Index for derivative should be 1, 2 or 3\n")
305 
306  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
307  Map<VectorXd> eigen_y(y.vector, y.vlen);
308 
309  Map<VectorXd> eigen_f(func.vector, func.vlen);
310 
311  SGVector<float64_t> r(func.vlen);
312  Map<VectorXd> eigen_r(r.vector, r.vlen);
313 
314  // compute s(f)=1./(1+exp(-f))
315  VectorXd eigen_s=(VectorXd::Ones(func.vlen)).cwiseQuotient((1.0+
316  (-eigen_f).array().exp()).matrix());
317 
318  // compute derivatives of log probability wrt f
319  if (i == 1)
320  {
321  // compute the first derivative: dlp=(y+1)/2-s(f)
322  eigen_r=(eigen_y.array()+1.0)/2.0-eigen_s.array();
323  }
324  else if (i == 2)
325  {
326  // compute the second derivative: d2lp=-s(f).*(1-s(f))
327  eigen_r=-eigen_s.array()*(1.0-eigen_s.array());
328  }
329  else if (i == 3)
330  {
331  // compute the third derivative: d2lp=-s(f).*(1-s(f)).*(1-2*s(f))
332  eigen_r=-eigen_s.array()*(1.0-eigen_s.array())*(1.0-2.0*eigen_s.array());
333  }
334  else
335  {
336  SG_ERROR("Invalid index for derivative\n")
337  }
338 
339  return r;
340 }
341 
343  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
344 {
346 
347  if (lab)
348  {
349  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
350  "Length of the vector of means (%d), length of the vector of "
351  "variances (%d) and number of labels (%d) should be the same\n",
352  mu.vlen, s2.vlen, lab->get_num_labels())
354  "Labels must be type of CBinaryLabels\n")
355 
356  y=((CBinaryLabels*)lab)->get_labels();
357  }
358  else
359  {
360  REQUIRE(mu.vlen==s2.vlen, "Length of the vector of means (%d) and "
361  "length of the vector of variances (%d) should be the same\n",
362  mu.vlen, s2.vlen)
363 
365  y.set_const(1.0);
366  }
367 
368  // create an object of normal pdf function
369  CNormalPDF* f=new CNormalPDF();
370 
371  // create an object of sigmoid function
372  CSigmoidFunction* g=new CSigmoidFunction();
373 
374  // create an object of product of sigmoid and normal pdf functions
375  CProductFunction* h=new CProductFunction(f, g);
376  SG_REF(h);
377 
378  // compute probabilities using numerical integration
380 
381  for (index_t i=0; i<mu.vlen; i++)
382  {
383  // set normal pdf parameters
384  f->set_mu(mu[i]);
385  f->set_sigma(CMath::sqrt(s2[i]));
386 
387  // set sigmoid parameters
388  g->set_a(y[i]);
389 
390  // evaluate integral on (-inf, inf)
393  }
394 
395  SG_UNREF(h);
396 
397  for (index_t i=0; i<r.vlen; i++)
398  r[i]=CMath::log(r[i]);
399 
400  return r;
401 }
402 
404  SGVector<float64_t> s2, const CLabels *lab, index_t i) const
405 {
406  // check the parameters
407  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
408  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
409  "Length of the vector of means (%d), length of the vector of "
410  "variances (%d) and number of labels (%d) should be the same\n",
411  mu.vlen, s2.vlen, lab->get_num_labels())
412  REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
414  "Labels must be type of CBinaryLabels\n")
415 
416  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
417 
418  // create an object of f(x)=N(x|mu,sigma^2)
419  CNormalPDF* f=new CNormalPDF(mu[i], CMath::sqrt(s2[i]));
420 
421  // create an object of g(x)=sigmoid(x)
422  CSigmoidFunction* g=new CSigmoidFunction(y[i]);
423 
424  // create an object of h(x)=N(x|mu,sigma^2)*sigmoid(x)
425  CProductFunction* h=new CProductFunction(f, g);
426 
427  // create an object of l(x)=x
428  CLinearFunction* l=new CLinearFunction();
429 
430  // create an object of k(x)=x*N(x|mu,sigma^2)*sigmoid(x)
431  CProductFunction* k=new CProductFunction(l, h);
432  SG_REF(k);
433 
434  // compute Z = \int N(x|mu,sigma)*sigmoid(a*x) dx
437 
438  // compute 1st moment: E[x] = Z^-1 * \int x*N(x|mu,sigma)*sigmoid(a*x)dx
441 
442  SG_UNREF(k);
443 
444  return Ex;
445 }
446 
448  SGVector<float64_t> s2, const CLabels *lab, index_t i) const
449 {
450  // check the parameters
451  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
452  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
453  "Length of the vector of means (%d), length of the vector of "
454  "variances (%d) and number of labels (%d) should be the same\n",
455  mu.vlen, s2.vlen, lab->get_num_labels())
456  REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
458  "Labels must be type of CBinaryLabels\n")
459 
460  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
461 
462  // create an object of f(x)=N(x|mu,sigma^2)
463  CNormalPDF* f=new CNormalPDF(mu[i], CMath::sqrt(s2[i]));
464 
465  // create an object of g(x)=sigmoid(a*x)
466  CSigmoidFunction* g=new CSigmoidFunction(y[i]);
467 
468  // create an object of h(x)=N(x|mu,sigma^2)*sigmoid(a*x)
469  CProductFunction* h=new CProductFunction(f, g);
470 
471  // create an object of l(x)=x
472  CLinearFunction* l=new CLinearFunction();
473 
474  // create an object of k(x)=x*N(x|mu,sigma^2)*sigmoid(a*x)
475  CProductFunction* k=new CProductFunction(l, h);
476  SG_REF(k);
477 
478  // create an object of q(x)=x^2
479  CQuadraticFunction* q=new CQuadraticFunction();
480 
481  // create an object of p(x)=x^2*N(x|mu,sigma^2)*sigmoid(x)
482  CProductFunction* p=new CProductFunction(q, h);
483  SG_REF(p);
484 
485  // compute Z = \int N(x|mu,sigma)*sigmoid(a*x) dx
488 
489  // compute 1st moment: E[x] = Z^-1 * \int x*N(x|mu,sigma)*sigmoid(a*x)dx
492 
493  // compute E[x^2] = Z^-1 * \int x^2*N(x|mu,sigma)*sigmoid(a*x)dx
496 
497  SG_UNREF(k);
498  SG_UNREF(p);
499 
500  // return 2nd moment: Var[x]=E[x^2]-E[x]^2
501  return Ex2-CMath::sq(Ex);;
502 }
503 }
504 
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:62
static float64_t integrate_quadgk(CFunction *f, float64_t a, float64_t b, float64_t abs_tol=1e-10, float64_t rel_tol=1e-5, uint32_t max_iter=1000, index_t sn=10)
The class Labels models labels, i.e. class assignments of objects.
Definition: Labels.h:43
static const float64_t INFTY
infinity
Definition: Math.h:2048
virtual int32_t get_num_labels() const =0
static T sq(T x)
Definition: Math.h:450
Definition: SGMatrix.h:20
#define SG_ERROR(...)
Definition: SGIO.h:129
#define REQUIRE(x,...)
Definition: SGIO.h:206
#define SG_REF(x)
Definition: SGObject.h:54
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
index_t vlen
Definition: SGVector.h:494
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
double float64_t
Definition: common.h:50
#define SG_UNREF(x)
Definition: SGObject.h:55
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:621
static float64_t log(float64_t v)
Definition: Math.h:922
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
static float32_t sqrt(float32_t x)
Definition: Math.h:459
The Likelihood model base class.
void set_const(T const_elem)
Definition: SGVector.cpp:150
static const float64_t PI
Definition: Math.h:2055

SHOGUN Machine Learning Toolbox - Documentation