SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LogitLikelihood.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2013 Roman Votyakov
8  */
9 
11 
12 #ifdef HAVE_EIGEN3
13 
18 
19 using namespace shogun;
20 using namespace Eigen;
21 
22 namespace shogun
23 {
24 
25 #ifndef DOXYGEN_SHOULD_SKIP_THIS
26 
28 class CNormalPDF : public CFunction
29 {
30 public:
32  CNormalPDF()
33  {
34  m_mu=0.0;
35  m_sigma=1.0;
36  }
37 
43  CNormalPDF(float64_t mu, float64_t sigma)
44  {
45  m_mu=mu;
46  m_sigma=sigma;
47  }
48 
53  void set_mu(float64_t mu) { m_mu=mu; }
54 
59  void set_sigma(float64_t sigma) { m_sigma=sigma; }
60 
67  virtual float64_t operator() (float64_t x)
68  {
69  return (1.0/(CMath::sqrt(2*CMath::PI)*m_sigma))*
70  CMath::exp(-CMath::sq(x-m_mu)/(2.0*CMath::sq(m_sigma)));
71  }
72 
73 private:
74  /* mean */
75  float64_t m_mu;
76 
77  /* standard deviation */
78  float64_t m_sigma;
79 };
80 
82 class CSigmoidFunction : public CFunction
83 {
84 public:
86  CSigmoidFunction()
87  {
88  m_a=0.0;
89  }
90 
95  CSigmoidFunction(float64_t a)
96  {
97  m_a=a;
98  }
99 
104  void set_a(float64_t a) { m_a=a; }
105 
112  virtual float64_t operator() (float64_t x)
113  {
114  return 1.0/(1.0+CMath::exp(-m_a*x));
115  }
116 
117 private:
119  float64_t m_a;
120 };
121 
123 class CProductFunction : public CFunction
124 {
125 public:
131  CProductFunction(CFunction* f, CFunction* g)
132  {
133  SG_REF(f);
134  SG_REF(g);
135  m_f=f;
136  m_g=g;
137  }
138 
139  virtual ~CProductFunction()
140  {
141  SG_UNREF(m_f);
142  SG_UNREF(m_g);
143  }
144 
151  virtual float64_t operator() (float64_t x)
152  {
153  return (*m_f)(x)*(*m_g)(x);
154  }
155 
156 private:
158  CFunction* m_f;
160  CFunction* m_g;
161 };
162 
164 class CLinearFunction : public CFunction
165 {
166 public:
168  CLinearFunction() { }
169 
170  virtual ~CLinearFunction() { }
171 
178  virtual float64_t operator() (float64_t x)
179  {
180  return x;
181  }
182 };
183 
185 class CQuadraticFunction : public CFunction
186 {
187 public:
189  CQuadraticFunction() { }
190 
191  virtual ~CQuadraticFunction() { }
192 
199  virtual float64_t operator() (float64_t x)
200  {
201  return CMath::sq(x);
202  }
203 };
204 
205 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
206 
208 {
209 }
210 
212 {
213 }
214 
216  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
217 {
219  Map<VectorXd> eigen_lp(lp.vector, lp.vlen);
220 
222  Map<VectorXd> eigen_r(r.vector, r.vlen);
223 
224  // evaluate predictive mean: ymu=2*p-1
225  // Note that the distribution is Bernoulli distribution with p(x=1)=p and
226  // p(x=-1)=(1-p)
227  // the mean of the Bernoulli distribution is 2*p-1
228  eigen_r=2.0*eigen_lp.array().exp()-1.0;
229 
230  return r;
231 }
232 
234  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
235 {
237  Map<VectorXd> eigen_lp(lp.vector, lp.vlen);
238 
240  Map<VectorXd> eigen_r(r.vector, r.vlen);
241 
242  // evaluate predictive variance: ys2=1-(2*p-1).^2
243  // Note that the distribution is Bernoulli distribution with p(x=1)=p and
244  // p(x=-1)=(1-p)
245  // the variance of the Bernoulli distribution is 1-(2*p-1).^2
246  eigen_r=1-(2.0*eigen_lp.array().exp()-1.0).square();
247 
248  return r;
249 }
250 
252  SGVector<float64_t> func) const
253 {
254  // check the parameters
255  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
257  "Labels must be type of CBinaryLabels\n")
258  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
259  "length of the function vector\n")
260 
261  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
262  Map<VectorXd> eigen_y(y.vector, y.vlen);
263 
264  Map<VectorXd> eigen_f(func.vector, func.vlen);
265 
266  SGVector<float64_t> r(func.vlen);
267  Map<VectorXd> eigen_r(r.vector, r.vlen);
268 
269  // compute log probability: -log(1+exp(-f.*y))
270  eigen_r=-(1.0+(-eigen_y.array()*eigen_f.array()).exp()).log();
271 
272  return r;
273 }
274 
276  const CLabels* lab, SGVector<float64_t> func, index_t i) const
277 {
278  // check the parameters
279  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
281  "Labels must be type of CBinaryLabels\n")
282  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
283  "length of the function vector\n")
284  REQUIRE(i>=1 && i<=3, "Index for derivative should be 1, 2 or 3\n")
285 
286  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
287  Map<VectorXd> eigen_y(y.vector, y.vlen);
288 
289  Map<VectorXd> eigen_f(func.vector, func.vlen);
290 
291  SGVector<float64_t> r(func.vlen);
292  Map<VectorXd> eigen_r(r.vector, r.vlen);
293 
294  // compute s(f)=1./(1+exp(-f))
295  VectorXd eigen_s=(VectorXd::Ones(func.vlen)).cwiseQuotient((1.0+
296  (-eigen_f).array().exp()).matrix());
297 
298  // compute derivatives of log probability wrt f
299  if (i == 1)
300  {
301  // compute the first derivative: dlp=(y+1)/2-s(f)
302  eigen_r=(eigen_y.array()+1.0)/2.0-eigen_s.array();
303  }
304  else if (i == 2)
305  {
306  // compute the second derivative: d2lp=-s(f).*(1-s(f))
307  eigen_r=-eigen_s.array()*(1.0-eigen_s.array());
308  }
309  else if (i == 3)
310  {
311  // compute the third derivative: d2lp=-s(f).*(1-s(f)).*(1-2*s(f))
312  eigen_r=-eigen_s.array()*(1.0-eigen_s.array())*(1.0-2.0*eigen_s.array());
313  }
314  else
315  {
316  SG_ERROR("Invalid index for derivative\n")
317  }
318 
319  return r;
320 }
321 
323  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
324 {
326 
327  if (lab)
328  {
329  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
330  "Length of the vector of means (%d), length of the vector of "
331  "variances (%d) and number of labels (%d) should be the same\n",
332  mu.vlen, s2.vlen, lab->get_num_labels())
334  "Labels must be type of CBinaryLabels\n")
335 
336  y=((CBinaryLabels*)lab)->get_labels();
337  }
338  else
339  {
340  REQUIRE(mu.vlen==s2.vlen, "Length of the vector of means (%d) and "
341  "length of the vector of variances (%d) should be the same\n",
342  mu.vlen, s2.vlen)
343 
345  y.set_const(1.0);
346  }
347 
348  // create an object of normal pdf function
349  CNormalPDF* f=new CNormalPDF();
350 
351  // create an object of sigmoid function
352  CSigmoidFunction* g=new CSigmoidFunction();
353 
354  // create an object of product of sigmoid and normal pdf functions
355  CProductFunction* h=new CProductFunction(f, g);
356  SG_REF(h);
357 
358  // compute probabilities using numerical integration
360 
361  for (index_t i=0; i<mu.vlen; i++)
362  {
363  // set normal pdf parameters
364  f->set_mu(mu[i]);
365  f->set_sigma(CMath::sqrt(s2[i]));
366 
367  // set sigmoid parameters
368  g->set_a(y[i]);
369 
370  // evaluate integral on (-inf, inf)
373  }
374 
375  SG_UNREF(h);
376 
377  r.log();
378 
379  return r;
380 }
381 
383  SGVector<float64_t> s2, const CLabels *lab, index_t i) const
384 {
385  // check the parameters
386  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
387  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
388  "Length of the vector of means (%d), length of the vector of "
389  "variances (%d) and number of labels (%d) should be the same\n",
390  mu.vlen, s2.vlen, lab->get_num_labels())
391  REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
393  "Labels must be type of CBinaryLabels\n")
394 
395  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
396 
397  // create an object of f(x)=N(x|mu,sigma^2)
398  CNormalPDF* f=new CNormalPDF(mu[i], CMath::sqrt(s2[i]));
399 
400  // create an object of g(x)=sigmoid(x)
401  CSigmoidFunction* g=new CSigmoidFunction(y[i]);
402 
403  // create an object of h(x)=N(x|mu,sigma^2)*sigmoid(x)
404  CProductFunction* h=new CProductFunction(f, g);
405 
406  // create an object of l(x)=x
407  CLinearFunction* l=new CLinearFunction();
408 
409  // create an object of k(x)=x*N(x|mu,sigma^2)*sigmoid(x)
410  CProductFunction* k=new CProductFunction(l, h);
411  SG_REF(k);
412 
413  // compute Z = \int N(x|mu,sigma)*sigmoid(a*x) dx
416 
417  // compute 1st moment: E[x] = Z^-1 * \int x*N(x|mu,sigma)*sigmoid(a*x)dx
420 
421  SG_UNREF(k);
422 
423  return Ex;
424 }
425 
427  SGVector<float64_t> s2, const CLabels *lab, index_t i) const
428 {
429  // check the parameters
430  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
431  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
432  "Length of the vector of means (%d), length of the vector of "
433  "variances (%d) and number of labels (%d) should be the same\n",
434  mu.vlen, s2.vlen, lab->get_num_labels())
435  REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
437  "Labels must be type of CBinaryLabels\n")
438 
439  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
440 
441  // create an object of f(x)=N(x|mu,sigma^2)
442  CNormalPDF* f=new CNormalPDF(mu[i], CMath::sqrt(s2[i]));
443 
444  // create an object of g(x)=sigmoid(a*x)
445  CSigmoidFunction* g=new CSigmoidFunction(y[i]);
446 
447  // create an object of h(x)=N(x|mu,sigma^2)*sigmoid(a*x)
448  CProductFunction* h=new CProductFunction(f, g);
449 
450  // create an object of l(x)=x
451  CLinearFunction* l=new CLinearFunction();
452 
453  // create an object of k(x)=x*N(x|mu,sigma^2)*sigmoid(a*x)
454  CProductFunction* k=new CProductFunction(l, h);
455  SG_REF(k);
456 
457  // create an object of q(x)=x^2
458  CQuadraticFunction* q=new CQuadraticFunction();
459 
460  // create an object of p(x)=x^2*N(x|mu,sigma^2)*sigmoid(x)
461  CProductFunction* p=new CProductFunction(q, h);
462  SG_REF(p);
463 
464  // compute Z = \int N(x|mu,sigma)*sigmoid(a*x) dx
467 
468  // compute 1st moment: E[x] = Z^-1 * \int x*N(x|mu,sigma)*sigmoid(a*x)dx
471 
472  // compute E[x^2] = Z^-1 * \int x^2*N(x|mu,sigma)*sigmoid(a*x)dx
475 
476  SG_UNREF(k);
477  SG_UNREF(p);
478 
479  // return 2nd moment: Var[x]=E[x^2]-E[x]^2
480  return Ex2-CMath::sq(Ex);;
481 }
482 }
483 
484 #endif /* HAVE_EIGEN3 */

SHOGUN Machine Learning Toolbox - Documentation