SHOGUN  6.1.3
StudentsTLikelihood.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (c) The Shogun Machine Learning Toolbox
3  * Written (W) 2013 Roman Votyakov
4  * Written (W) 2012 Jacob Walker
5  * All rights reserved.
6  *
7  * Redistribution and use in source and binary forms, with or without
8  * modification, are permitted provided that the following conditions are met:
9  *
10  * 1. Redistributions of source code must retain the above copyright notice, this
11  * list of conditions and the following disclaimer.
12  * 2. Redistributions in binary form must reproduce the above copyright notice,
13  * this list of conditions and the following disclaimer in the documentation
14  * and/or other materials provided with the distribution.
15  *
16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
18  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
20  * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26  *
27  * The views and conclusions contained in the software and documentation are those
28  * of the authors and should not be interpreted as representing official policies,
29  * either expressed or implied, of the Shogun Development Team.
30  *
31  * Adapted from the GPML toolbox, specifically likT.m
32  * http://www.gaussianprocess.org/gpml/code/matlab/doc/
33  */
35 
36 
38 #ifdef USE_GPL_SHOGUN
39 #include <shogun/mathematics/Integration.h>
40 #endif //USE_GPL_SHOGUN
45 
46 using namespace shogun;
47 using namespace Eigen;
48 
49 namespace shogun
50 {
51 
52 #ifndef DOXYGEN_SHOULD_SKIP_THIS
53 
55 class CNormalPDF : public CFunction
56 {
57 public:
59  CNormalPDF()
60  {
61  m_mu=0.0;
62  m_sigma=1.0;
63  }
64 
70  CNormalPDF(float64_t mu, float64_t sigma)
71  {
72  m_mu=mu;
73  m_sigma=sigma;
74  }
75 
80  void set_mu(float64_t mu) { m_mu=mu; }
81 
86  void set_sigma(float64_t sigma) { m_sigma=sigma; }
87 
94  virtual float64_t operator() (float64_t x)
95  {
96  return (1.0/(CMath::sqrt(2*CMath::PI)*m_sigma))*
97  CMath::exp(-CMath::sq(x-m_mu)/(2.0*CMath::sq(m_sigma)));
98  }
99 
100 private:
101  /* mean */
102  float64_t m_mu;
103 
104  /* standard deviation */
105  float64_t m_sigma;
106 };
107 
111 class CStudentsTPDF : public CFunction
112 {
113 public:
115  CStudentsTPDF()
116  {
117  m_sigma=1.0;
118  m_mu=0.0;
119  m_nu=3.0;
120  }
121 
128  CStudentsTPDF(float64_t sigma, float64_t nu, float64_t mu)
129  {
130  m_sigma=sigma;
131  m_mu=mu;
132  m_nu=nu;
133  }
134 
139  void set_sigma(float64_t sigma) { m_sigma=sigma; }
140 
145  void set_mu(float64_t mu) { m_mu=mu; }
146 
151  void set_nu(float64_t nu) { m_nu=nu; }
152 
160  virtual float64_t operator() (float64_t x)
161  {
162  float64_t lZ=CStatistics::lgamma((m_nu+1.0)/2.0)-CStatistics::lgamma(m_nu/2.0)-
163  CMath::log(m_nu*CMath::PI*CMath::sq(m_sigma))/2.0;
164  return CMath::exp(lZ-((m_nu+1.0)/2.0)*CMath::log(1.0+CMath::sq(x-m_mu)/
165  (m_nu*CMath::sq(m_sigma))));
166  }
167 
168 private:
170  float64_t m_sigma;
171 
173  float64_t m_mu;
174 
176  float64_t m_nu;
177 };
178 
180 class CProductFunction : public CFunction
181 {
182 public:
188  CProductFunction(CFunction* f, CFunction* g)
189  {
190  SG_REF(f);
191  SG_REF(g);
192  m_f=f;
193  m_g=g;
194  }
195 
196  virtual ~CProductFunction()
197  {
198  SG_UNREF(m_f);
199  SG_UNREF(m_g);
200  }
201 
208  virtual float64_t operator() (float64_t x)
209  {
210  return (*m_f)(x)*(*m_g)(x);
211  }
212 
213 private:
215  CFunction* m_f;
217  CFunction* m_g;
218 };
219 
221 class CLinearFunction : public CFunction
222 {
223 public:
225  CLinearFunction() { }
226 
227  virtual ~CLinearFunction() { }
228 
235  virtual float64_t operator() (float64_t x)
236  {
237  return x;
238  }
239 };
240 
242 class CQuadraticFunction : public CFunction
243 {
244 public:
246  CQuadraticFunction() { }
247 
248  virtual ~CQuadraticFunction() { }
249 
256  virtual float64_t operator() (float64_t x)
257  {
258  return CMath::sq(x);
259  }
260 };
261 
262 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
263 
265 {
266  init();
267 }
268 
270  : CLikelihoodModel()
271 {
272  init();
273  set_sigma(sigma);
275 }
276 
277 void CStudentsTLikelihood::init()
278 {
279  m_log_sigma=0.0;
280  m_log_df=CMath::log(2.0);
281  SG_ADD(&m_log_df, "log_df", "Degrees of freedom in log domain", MS_AVAILABLE, GRADIENT_AVAILABLE);
282  SG_ADD(&m_log_sigma, "log_sigma", "Scale parameter in log domain", MS_AVAILABLE, GRADIENT_AVAILABLE);
283 }
284 
286 {
287 }
288 
290  CLikelihoodModel* lik)
291 {
292  ASSERT(lik!=NULL);
293 
294  if (lik->get_model_type()!=LT_STUDENTST)
295  SG_SERROR("Provided likelihood is not of type CStudentsTLikelihood!\n")
296 
297  SG_REF(lik);
298  return (CStudentsTLikelihood*)lik;
299 }
300 
302  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
303 {
304  return SGVector<float64_t>(mu);
305 }
306 
308  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
309 {
310  SGVector<float64_t> result(s2);
311  Map<VectorXd> eigen_result(result.vector, result.vlen);
313  if (df<2.0)
314  eigen_result=CMath::INFTY*VectorXd::Ones(result.vlen);
315  else
316  {
317  eigen_result+=CMath::exp(m_log_sigma*2.0)*df/(df-2.0)*
318  VectorXd::Ones(result.vlen);
319  }
320 
321  return result;
322 }
323 
325  SGVector<float64_t> func) const
326 {
327  // check the parameters
328  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
330  "Labels must be type of CRegressionLabels\n")
331  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
332  "length of the function vector\n")
333 
334  Map<VectorXd> eigen_f(func.vector, func.vlen);
335 
336  SGVector<float64_t> r(func.vlen);
337  Map<VectorXd> eigen_r(r.vector, r.vlen);
338 
339  SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels();
340  Map<VectorXd> eigen_y(y.vector, y.vlen);
341 
343  // compute lZ=log(gamma(df/2+1/2))-log(gamma(df/2))-log(df*pi*sigma^2)/2
344  VectorXd eigen_lZ=(CStatistics::lgamma(df/2.0+0.5)-
345  CStatistics::lgamma(df/2.0)-log(df*CMath::PI*CMath::exp(m_log_sigma*2.0))/2.0)*
346  VectorXd::Ones(r.vlen);
347 
348  // compute log probability: lp=lZ-(df+1)*log(1+(y-f).^2./(df*sigma^2))/2
349  eigen_r=eigen_y-eigen_f;
350  eigen_r=eigen_r.cwiseProduct(eigen_r)/(df*CMath::exp(m_log_sigma*2.0));
351  eigen_r=eigen_lZ-(df+1)*
352  (eigen_r+VectorXd::Ones(r.vlen)).array().log().matrix()/2.0;
353 
354  return r;
355 }
356 
358  const CLabels* lab, SGVector<float64_t> func, index_t i) const
359 {
360  // check the parameters
361  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
363  "Labels must be type of CRegressionLabels\n")
364  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
365  "length of the function vector\n")
366  REQUIRE(i>=1 && i<=3, "Index for derivative should be 1, 2 or 3\n")
367 
368  Map<VectorXd> eigen_f(func.vector, func.vlen);
369 
370  SGVector<float64_t> r(func.vlen);
371  Map<VectorXd> eigen_r(r.vector, r.vlen);
372 
373  SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels();
374  Map<VectorXd> eigen_y(y.vector, y.vlen);
375 
376  // compute r=y-f, r2=r.^2
377  eigen_r=eigen_y-eigen_f;
378  VectorXd eigen_r2=eigen_r.cwiseProduct(eigen_r);
380 
381  // compute a=(y-f).^2+df*sigma^2
382  VectorXd a=eigen_r2+VectorXd::Ones(r.vlen)*df*CMath::exp(m_log_sigma*2.0);
383 
384  if (i==1)
385  {
386  // compute first derivative of log probability wrt f:
387  // dlp=(df+1)*(y-f)./a
388  eigen_r=(df+1)*eigen_r.cwiseQuotient(a);
389  }
390  else if (i==2)
391  {
392  // compute second derivative of log probability wrt f:
393  // d2lp=(df+1)*((y-f)^2-df*sigma^2)./a.^2
394  VectorXd b=eigen_r2-VectorXd::Ones(r.vlen)*df*CMath::exp(m_log_sigma*2.0);
395 
396  eigen_r=(df+1)*b.cwiseQuotient(a.cwiseProduct(a));
397  }
398  else if (i==3)
399  {
400  // compute third derivative of log probability wrt f:
401  // d3lp=(f+1)*2*(y-f).*((y-f)^2-3*f*sigma^2)./a.^3
402  VectorXd c=eigen_r2-VectorXd::Ones(r.vlen)*3*df*CMath::exp(m_log_sigma*2.0);
403  VectorXd a2=a.cwiseProduct(a);
404 
405  eigen_r=(df+1)*2*eigen_r.cwiseProduct(c).cwiseQuotient(
406  a2.cwiseProduct(a));
407  }
408 
409  return r;
410 }
411 
413  SGVector<float64_t> func, const TParameter* param) const
414 {
415  // check the parameters
416  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
418  "Labels must be type of CRegressionLabels\n")
419  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
420  "length of the function vector\n")
421 
422  SGVector<float64_t> r(func.vlen);
423  Map<VectorXd> eigen_r(r.vector, r.vlen);
424 
425  Map<VectorXd> eigen_f(func.vector, func.vlen);
426 
427  SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels();
428  Map<VectorXd> eigen_y(y.vector, y.vlen);
429 
430  // compute r=y-f and r2=(y-f).^2
431  eigen_r=eigen_y-eigen_f;
432  VectorXd eigen_r2=eigen_r.cwiseProduct(eigen_r);
434 
435  if (!strcmp(param->m_name, "log_df"))
436  {
437  // compute derivative of log probability wrt df:
438  // lp_ddf=df*(dloggamma(df/2+1/2)-dloggamma(df/2))/2-1/2-
439  // df*log(1+r2/(df*sigma^2))/2 +(df/2+1/2)*r2./(df*sigma^2+r2)
440  eigen_r=(df*(CStatistics::dlgamma(df*0.5+0.5)-
441  CStatistics::dlgamma(df*0.5))*0.5-0.5)*VectorXd::Ones(r.vlen);
442 
443  eigen_r-=df*(VectorXd::Ones(r.vlen)+
444  eigen_r2/(df*CMath::exp(m_log_sigma*2.0))).array().log().matrix()/2.0;
445 
446  eigen_r+=(df/2.0+0.5)*eigen_r2.cwiseQuotient(
447  eigen_r2+VectorXd::Ones(r.vlen)*(df*CMath::exp(m_log_sigma*2.0)));
448 
449  eigen_r*=(1.0-1.0/df);
450 
451  return r;
452  }
453  else if (!strcmp(param->m_name, "log_sigma"))
454  {
455  // compute derivative of log probability wrt sigma:
456  // lp_dsigma=(df+1)*r2./a-1
457  eigen_r=(df+1)*eigen_r2.cwiseQuotient(eigen_r2+
458  VectorXd::Ones(r.vlen)*(df*CMath::exp(m_log_sigma*2.0)));
459  eigen_r-=VectorXd::Ones(r.vlen);
460 
461  return r;
462  }
463 
464  return SGVector<float64_t>();
465 }
466 
468  SGVector<float64_t> func, const TParameter* param) const
469 {
470  // check the parameters
471  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
473  "Labels must be type of CRegressionLabels\n")
474  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
475  "length of the function vector\n")
476 
477  SGVector<float64_t> r(func.vlen);
478  Map<VectorXd> eigen_r(r.vector, r.vlen);
479 
480  Map<VectorXd> eigen_f(func.vector, func.vlen);
481 
482  SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels();
483  Map<VectorXd> eigen_y(y.vector, y.vlen);
484 
485  // compute r=y-f and r2=(y-f).^2
486  eigen_r=eigen_y-eigen_f;
487  VectorXd eigen_r2=eigen_r.cwiseProduct(eigen_r);
489 
490  // compute a=r+sigma^2*df and a2=a.^2
491  VectorXd a=eigen_r2+CMath::exp(m_log_sigma*2.0)*df*VectorXd::Ones(r.vlen);
492  VectorXd a2=a.cwiseProduct(a);
493 
494  if (!strcmp(param->m_name, "log_df"))
495  {
496  // compute derivative of first derivative of log probability wrt df:
497  // dlp_ddf=df*r.*(a-sigma^2*(df+1))./a2
498  eigen_r=df*eigen_r.cwiseProduct(a-CMath::exp(m_log_sigma*2.0)*(df+1.0)*
499  VectorXd::Ones(r.vlen)).cwiseQuotient(a2);
500  eigen_r*=(1.0-1.0/df);
501 
502  return r;
503  }
504  else if (!strcmp(param->m_name, "log_sigma"))
505  {
506  // compute derivative of first derivative of log probability wrt sigma:
507  // dlp_dsigma=-(df+1)*2*df*sigma^2*r./a2
508  eigen_r=-(df+1.0)*2*df*CMath::exp(m_log_sigma*2.0)*
509  eigen_r.cwiseQuotient(a2);
510 
511  return r;
512  }
513 
514  return SGVector<float64_t>();
515 }
516 
518  SGVector<float64_t> func, const TParameter* param) const
519 {
520  // check the parameters
521  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
523  "Labels must be type of CRegressionLabels\n")
524  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
525  "length of the function vector\n")
526 
527  SGVector<float64_t> r(func.vlen);
528  Map<VectorXd> eigen_r(r.vector, r.vlen);
529 
530  Map<VectorXd> eigen_f(func.vector, func.vlen);
531 
532  SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels();
533  Map<VectorXd> eigen_y(y.vector, y.vlen);
534 
535  // compute r=y-f and r2=(y-f).^2
536  eigen_r=eigen_y-eigen_f;
537  VectorXd eigen_r2=eigen_r.cwiseProduct(eigen_r);
539 
540  // compute a=r+sigma^2*df and a3=a.^3
541  VectorXd a=eigen_r2+CMath::exp(m_log_sigma*2.0)*df*VectorXd::Ones(r.vlen);
542  VectorXd a3=(a.cwiseProduct(a)).cwiseProduct(a);
543 
544  if (!strcmp(param->m_name, "log_df"))
545  {
546  // compute derivative of second derivative of log probability wrt df:
547  // d2lp_ddf=df*(r2.*(r2-3*sigma^2*(1+df))+df*sigma^4)./a3
548  float64_t sigma2=CMath::exp(m_log_sigma*2.0);
549 
550  eigen_r=df*(eigen_r2.cwiseProduct(eigen_r2-3*sigma2*(1.0+df)*
551  VectorXd::Ones(r.vlen))+(df*CMath::sq(sigma2))*VectorXd::Ones(r.vlen));
552  eigen_r=eigen_r.cwiseQuotient(a3);
553 
554  eigen_r*=(1.0-1.0/df);
555 
556  return r;
557  }
558  else if (!strcmp(param->m_name, "log_sigma"))
559  {
560  // compute derivative of second derivative of log probability wrt sigma:
561  // d2lp_dsigma=(df+1)*2*df*sigma^2*(a-4*r2)./a3
562  eigen_r=(df+1.0)*2*df*CMath::exp(m_log_sigma*2.0)*
563  (a-4.0*eigen_r2).cwiseQuotient(a3);
564 
565  return r;
566  }
567 
568  return SGVector<float64_t>();
569 }
570 
572  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
573 {
575 
576  if (lab)
577  {
578  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
579  "Length of the vector of means (%d), length of the vector of "
580  "variances (%d) and number of labels (%d) should be the same\n",
581  mu.vlen, s2.vlen, lab->get_num_labels())
583  "Labels must be type of CRegressionLabels\n")
584 
585  y=((CRegressionLabels*)lab)->get_labels();
586  }
587  else
588  {
589  REQUIRE(mu.vlen==s2.vlen, "Length of the vector of means (%d) and "
590  "length of the vector of variances (%d) should be the same\n",
591  mu.vlen, s2.vlen)
592 
594  y.set_const(1.0);
595  }
596 
597  // create an object of normal pdf
598  CNormalPDF* f=new CNormalPDF();
599 
600  // create an object of Student's t pdf
601  CStudentsTPDF* g=new CStudentsTPDF();
602 
603  g->set_nu(get_degrees_freedom());
604  g->set_sigma(CMath::exp(m_log_sigma));
605 
606  // create an object of product of Student's-t pdf and normal pdf
607  CProductFunction* h=new CProductFunction(f, g);
608  SG_REF(h);
609 
610  // compute probabilities using numerical integration
612 
613  for (index_t i=0; i<mu.vlen; i++)
614  {
615  // set normal pdf parameters
616  f->set_mu(mu[i]);
617  f->set_sigma(CMath::sqrt(s2[i]));
618 
619  // set Stundent's-t pdf parameters
620  g->set_mu(y[i]);
621 
622 #ifdef USE_GPL_SHOGUN
623  // evaluate integral on (-inf, inf)
624  r[i]=CIntegration::integrate_quadgk(h, -CMath::INFTY, mu[i])+
625  CIntegration::integrate_quadgk(h, mu[i], CMath::INFTY);
626 #else
627  SG_ERROR("StudentsT likelihood moments only supported under GPL.\n")
628 #endif //USE_GPL_SHOGUN
629  }
630 
631  SG_UNREF(h);
632 
633  for (index_t i=0; i<r.vlen; i++)
634  r[i]=CMath::log(r[i]);
635 
636  return r;
637 }
638 
640  SGVector<float64_t> s2, const CLabels *lab, index_t i) const
641 {
642  // check the parameters
643  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
644  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
645  "Length of the vector of means (%d), length of the vector of "
646  "variances (%d) and number of labels (%d) should be the same\n",
647  mu.vlen, s2.vlen, lab->get_num_labels())
648  REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
650  "Labels must be type of CRegressionLabels\n")
651 
652  SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels();
653 
654  // create an object of normal pdf
655  CNormalPDF* f=new CNormalPDF(mu[i], CMath::sqrt(s2[i]));
656 
657  // create an object of Student's t pdf
658  CStudentsTPDF* g=new CStudentsTPDF(CMath::exp(m_log_sigma), get_degrees_freedom(), y[i]);
659 
660  // create an object of h(x)=N(x|mu,sigma)*t(x|mu,sigma,nu)
661  CProductFunction* h=new CProductFunction(f, g);
662 
663  // create an object of k(x)=x*N(x|mu,sigma)*t(x|mu,sigma,nu)
664  CProductFunction* k=new CProductFunction(new CLinearFunction(), h);
665  SG_REF(k);
666 
667  float64_t Ex=0;
668 #ifdef USE_GPL_SHOGUN
669  // compute Z = \int N(x|mu,sigma)*t(x|mu,sigma,nu) dx
670  float64_t Z=CIntegration::integrate_quadgk(h, -CMath::INFTY, mu[i])+
671  CIntegration::integrate_quadgk(h, mu[i], CMath::INFTY);
672 
673  // compute 1st moment:
674  // E[x] = Z^-1 * \int x*N(x|mu,sigma)*t(x|mu,sigma,nu)dx
675  Ex=(CIntegration::integrate_quadgk(k, -CMath::INFTY, mu[i])+
676  CIntegration::integrate_quadgk(k, mu[i], CMath::INFTY))/Z;
677 #else
678  SG_ERROR("StudentsT likelihood moments only supported under GPL.\n")
679 #endif //USE_GPL_SHOGUN
680  SG_UNREF(k);
681 
682  return Ex;
683 }
684 
686  SGVector<float64_t> s2, const CLabels *lab, index_t i) const
687 {
688  // check the parameters
689  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
690  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
691  "Length of the vector of means (%d), length of the vector of "
692  "variances (%d) and number of labels (%d) should be the same\n",
693  mu.vlen, s2.vlen, lab->get_num_labels())
694  REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
696  "Labels must be type of CRegressionLabels\n")
697 
698  SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels();
699 
700  // create an object of normal pdf
701  CNormalPDF* f=new CNormalPDF(mu[i], CMath::sqrt(s2[i]));
702 
703  // create an object of Student's t pdf
704  CStudentsTPDF* g=new CStudentsTPDF(CMath::exp(m_log_sigma), get_degrees_freedom(), y[i]);
705 
706  // create an object of h(x)=N(x|mu,sigma)*t(x|mu,sigma,nu)
707  CProductFunction* h=new CProductFunction(f, g);
708 
709  // create an object of k(x)=x*N(x|mu,sigma)*t(x|mu,sigma,nu)
710  CProductFunction* k=new CProductFunction(new CLinearFunction(), h);
711  SG_REF(k);
712 
713  // create an object of p(x)=x^2*N(x|mu,sigma^2)*t(x|mu,sigma,nu)
714  CProductFunction* p=new CProductFunction(new CQuadraticFunction(), h);
715  SG_REF(p);
716 
717  float64_t Ex=0;
718  float64_t Ex2=0;
719 #ifdef USE_GPL_SHOGUN
720  // compute Z = \int N(x|mu,sigma)*t(x|mu,sigma,nu) dx
721  float64_t Z=CIntegration::integrate_quadgk(h, -CMath::INFTY, mu[i])+
722  CIntegration::integrate_quadgk(h, mu[i], CMath::INFTY);
723 
724  // compute 1st moment:
725  // E[x] = Z^-1 * \int x*N(x|mu,sigma)*t(x|mu,sigma,nu)dx
726  Ex=(CIntegration::integrate_quadgk(k, -CMath::INFTY, mu[i])+
727  CIntegration::integrate_quadgk(k, mu[i], CMath::INFTY))/Z;
728 
729  // compute E[x^2] = Z^-1 * \int x^2*N(x|mu,sigma)*t(x|mu,sigma,nu)dx
730  Ex2=(CIntegration::integrate_quadgk(p, -CMath::INFTY, mu[i])+
731  CIntegration::integrate_quadgk(p, mu[i], CMath::INFTY))/Z;
732 #else
734 #endif //USE_GPL_SHOGUN
735  SG_UNREF(k);
736  SG_UNREF(p);
737 
738  // return 2nd moment: Var[x]=E[x^2]-E[x]^2
739  return Ex2-CMath::sq(Ex);
740 }
741 }
742 
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
Real Labels are real-valued labels.
virtual float64_t get_second_moment(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab, index_t i) const
static float64_t lgamma(float64_t x)
Definition: Statistics.h:131
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 SGVector< float64_t > get_second_derivative(const CLabels *lab, SGVector< float64_t > func, const TParameter *param) const
virtual int32_t get_num_labels() const =0
real valued labels (e.g. for regression, classifier outputs)
Definition: LabelTypes.h:22
static T sqrt(T x)
Definition: Math.h:428
virtual SGVector< float64_t > get_first_derivative(const CLabels *lab, SGVector< float64_t > func, const TParameter *param) const
static T sq(T x)
Definition: Math.h:418
Definition: SGMatrix.h:25
virtual ELikelihoodModelType get_model_type() const
parameter struct
#define SG_ERROR(...)
Definition: SGIO.h:128
#define REQUIRE(x,...)
Definition: SGIO.h:181
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_log_probability_f(const CLabels *lab, SGVector< float64_t > func) const
float64_t get_degrees_freedom() const
virtual float64_t get_first_moment(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab, index_t i) const
virtual SGVector< float64_t > get_log_probability_derivative_f(const CLabels *lab, SGVector< float64_t > func, index_t i) const
#define SG_REF(x)
Definition: SGObject.h:52
Class of a function of one variable.
Definition: Function.h:22
#define ASSERT(x)
Definition: SGIO.h:176
#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
void set_sigma(float64_t sigma)
Class that models a Student&#39;s-t likelihood.
#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 dlgamma(float64_t x)
Definition: Statistics.cpp:561
void set_degrees_freedom(float64_t df)
#define SG_SERROR(...)
Definition: SGIO.h:164
static float64_t exp(float64_t x)
Definition: Math.h:551
static float64_t log(float64_t v)
Definition: Math.h:714
#define SG_ADD(...)
Definition: SGObject.h:93
static CStudentsTLikelihood * obtain_from_generic(CLikelihoodModel *likelihood)
virtual SGVector< float64_t > get_third_derivative(const CLabels *lab, SGVector< float64_t > func, const TParameter *param) const
The Likelihood model base class.
virtual SGVector< float64_t > get_log_zeroth_moments(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab) const
index_t vlen
Definition: SGVector.h:571
static const float64_t PI
Definition: Math.h:1875

SHOGUN Machine Learning Toolbox - Documentation