SHOGUN  4.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
GaussianLikelihood.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  */
33 
34 #ifdef HAVE_EIGEN3
35 
38 
39 using namespace shogun;
40 using namespace Eigen;
41 
43 {
44  init();
45 }
46 
48 {
49  init();
50  set_sigma(sigma);
51 }
52 
53 void CGaussianLikelihood::init()
54 {
55  m_log_sigma=0.0;
56  SG_ADD(&m_log_sigma, "log_sigma", "Observation noise in log domain", MS_AVAILABLE, GRADIENT_AVAILABLE);
57 }
58 
60 {
61 }
62 
64  CLikelihoodModel* lik)
65 {
66  ASSERT(lik!=NULL);
67 
68  if (lik->get_model_type()!=LT_GAUSSIAN)
69  SG_SERROR("Provided likelihood is not of type CGaussianLikelihood!\n")
70 
71  SG_REF(lik);
72  return (CGaussianLikelihood*)lik;
73 }
74 
76  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
77 {
78  return SGVector<float64_t>(mu);
79 }
80 
82  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
83 {
84  SGVector<float64_t> result(s2);
85  Map<VectorXd> eigen_result(result.vector, result.vlen);
86 
87  eigen_result=eigen_result.array()+CMath::exp(m_log_sigma*2.0);
88 
89  return result;
90 }
91 
93  SGVector<float64_t> func) const
94 {
95  // check the parameters
96  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
98  "Labels must be type of CRegressionLabels\n")
99  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
100  "length of the function vector\n")
101 
102  Map<VectorXd> eigen_f(func.vector, func.vlen);
103 
104  SGVector<float64_t> result(func.vlen);
105  Map<VectorXd> eigen_result(result.vector, result.vlen);
106 
107  SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels();
108  Map<VectorXd> eigen_y(y.vector, y.vlen);
109 
110  // compute log probability: lp=-(y-f).^2./sigma^2/2-log(2*pi*sigma^2)/2
111  eigen_result=eigen_y-eigen_f;
112  eigen_result=-eigen_result.cwiseProduct(eigen_result)/(2.0*CMath::exp(m_log_sigma*2.0))-
113  VectorXd::Ones(result.vlen)*log(2.0*CMath::PI*CMath::exp(m_log_sigma*2.0))/2.0;
114 
115  return result;
116 }
117 
119  const CLabels* lab, SGVector<float64_t> func, index_t i) const
120 {
121  // check the parameters
122  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
124  "Labels must be type of CRegressionLabels\n")
125  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
126  "length of the function vector\n")
127  REQUIRE(i>=1 && i<=3, "Index for derivative should be 1, 2 or 3\n")
128 
129  Map<VectorXd> eigen_f(func.vector, func.vlen);
130 
131  SGVector<float64_t> result(func.vlen);
132  Map<VectorXd> eigen_result(result.vector, result.vlen);
133 
134  SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels();
135  Map<VectorXd> eigen_y(y.vector, y.vlen);
136 
137  // set result=y-f
138  eigen_result=eigen_y-eigen_f;
139 
140  // compute derivatives of log probability wrt f
141  if (i == 1)
142  eigen_result/=CMath::exp(m_log_sigma*2.0);
143  else if (i == 2)
144  eigen_result=-VectorXd::Ones(result.vlen)/CMath::exp(m_log_sigma*2.0);
145  else if (i == 3)
146  eigen_result=VectorXd::Zero(result.vlen);
147 
148  return result;
149 }
150 
152  SGVector<float64_t> func, const TParameter* param) const
153 {
154  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
156  "Labels must be type of CRegressionLabels\n")
157  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
158  "length of the function vector\n")
159 
160  Map<VectorXd> eigen_f(func.vector, func.vlen);
161 
162  SGVector<float64_t> result(func.vlen);
163  Map<VectorXd> eigen_result(result.vector, result.vlen);
164 
165  if (strcmp(param->m_name, "log_sigma"))
166  return SGVector<float64_t>();
167 
168  SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels();
169  Map<VectorXd> eigen_y(y.vector, y.vlen);
170 
171  // compute derivative of log probability wrt log_sigma:
172  // dlp_dlogsigma
173  // lp_dsigma=(y-f).^2/sigma^2-1
174  eigen_result=eigen_y-eigen_f;
175  eigen_result=eigen_result.cwiseProduct(eigen_result)/CMath::exp(m_log_sigma*2.0);
176  eigen_result-=VectorXd::Ones(result.vlen);
177 
178  return result;
179 }
180 
182  SGVector<float64_t> func, const TParameter* param) const
183 {
184  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
186  "Labels must be type of CRegressionLabels\n")
187  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
188  "length of the function vector\n")
189 
190  if (strcmp(param->m_name, "log_sigma"))
191  return SGVector<float64_t>();
192 
193  Map<VectorXd> eigen_f(func.vector, func.vlen);
194 
195  SGVector<float64_t> result(func.vlen);
196  Map<VectorXd> eigen_result(result.vector, result.vlen);
197 
198  SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels();
199  Map<VectorXd> eigen_y(y.vector, y.vlen);
200 
201  // compute derivative of (the first log_sigma derivative of log probability) wrt f:
202  // d2lp_dlogsigma_df == d2lp_df_dlogsigma
203  // dlp_dsigma=2*(f-y)/sigma^2
204  eigen_result=2.0*(eigen_f-eigen_y)/CMath::exp(m_log_sigma*2.0);
205 
206  return result;
207 }
208 
210  SGVector<float64_t> func, const TParameter* param) const
211 {
212  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
214  "Labels must be type of CRegressionLabels\n")
215  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
216  "length of the function vector\n")
217 
218  if (strcmp(param->m_name, "log_sigma"))
219  return SGVector<float64_t>();
220 
221  Map<VectorXd> eigen_f(func.vector, func.vlen);
222 
223  SGVector<float64_t> result(func.vlen);
224  Map<VectorXd> eigen_result(result.vector, result.vlen);
225 
226  // compute derivative of (the derivative of the first log_sigma derivative of log probability) wrt f:
227  // d3lp_dlogsigma_df_df == d3lp_df_df_dlogsigma
228  // d2lp_dsigma=2/sigma^2
229  eigen_result=2.0*VectorXd::Ones(result.vlen)/CMath::exp(m_log_sigma*2.0);
230 
231  return result;
232 }
233 
235  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels *lab) const
236 {
238 
239  if (lab)
240  {
241  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
242  "Length of the vector of means (%d), length of the vector of "
243  "variances (%d) and number of labels (%d) should be the same\n",
244  mu.vlen, s2.vlen, lab->get_num_labels())
246  "Labels must be type of CRegressionLabels\n")
247 
248  y=((CRegressionLabels*)lab)->get_labels();
249  }
250  else
251  {
252  REQUIRE(mu.vlen==s2.vlen, "Length of the vector of means (%d) and "
253  "length of the vector of variances (%d) should be the same\n",
254  mu.vlen, s2.vlen)
255 
257  y.set_const(1.0);
258  }
259 
260  // create eigen representation of y, mu and s2
261  Map<VectorXd> eigen_mu(mu.vector, mu.vlen);
262  Map<VectorXd> eigen_s2(s2.vector, s2.vlen);
263  Map<VectorXd> eigen_y(y.vector, y.vlen);
264 
265  SGVector<float64_t> result(mu.vlen);
266  Map<VectorXd> eigen_result(result.vector, result.vlen);
267 
268  // compule lZ=-(y-mu).^2./(sn2+s2)/2-log(2*pi*(sn2+s2))/2
269  eigen_s2=eigen_s2.array()+CMath::exp(m_log_sigma*2.0);
270  eigen_result=-(eigen_y-eigen_mu).array().square()/(2.0*eigen_s2.array())-
271  (2.0*CMath::PI*eigen_s2.array()).log()/2.0;
272 
273  return result;
274 }
275 
277  SGVector<float64_t> s2, const CLabels *lab, index_t i) const
278 {
279  // check the parameters
280  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
281  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
282  "Length of the vector of means (%d), length of the vector of "
283  "variances (%d) and number of labels (%d) should be the same\n",
284  mu.vlen, s2.vlen, lab->get_num_labels())
285  REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
287  "Labels must be type of CRegressionLabels\n")
288 
289  SGVector<float64_t> y=((CRegressionLabels*)lab)->get_labels();
290 
291  // compute 1st moment
292  float64_t Ex=mu[i]+s2[i]*(y[i]-mu[i])/(CMath::exp(m_log_sigma*2.0)+s2[i]);
293 
294  return Ex;
295 }
296 
298  SGVector<float64_t> s2, const CLabels *lab, index_t i) const
299 {
300  // check the parameters
301  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
302  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
303  "Length of the vector of means (%d), length of the vector of "
304  "variances (%d) and number of labels (%d) should be the same\n",
305  mu.vlen, s2.vlen, lab->get_num_labels())
306  REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
308  "Labels must be type of CRegressionLabels\n")
309 
310  // compute 2nd moment
311  float64_t Var=s2[i]-CMath::sq(s2[i])/(CMath::exp(m_log_sigma*2.0)+s2[i]);
312 
313  return Var;
314 }
315 
316 #endif /* HAVE_EIGEN3 */
virtual ELabelType get_label_type() const =0
Class that models Gaussian likelihood.
Real Labels are real-valued labels.
virtual SGVector< float64_t > get_log_probability_f(const CLabels *lab, SGVector< float64_t > func) const
int32_t index_t
Definition: common.h:62
The class Labels models labels, i.e. class assignments of objects.
Definition: Labels.h:43
virtual int32_t get_num_labels() const =0
real valued labels (e.g. for regression, classifier outputs)
Definition: LabelTypes.h:22
static T sq(T x)
Definition: Math.h:450
Definition: SGMatrix.h:20
virtual ELikelihoodModelType get_model_type() const
parameter struct
#define REQUIRE(x,...)
Definition: SGIO.h:206
void set_sigma(float64_t sigma)
#define SG_REF(x)
Definition: SGObject.h:51
virtual float64_t get_second_moment(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab, index_t i) const
virtual SGVector< float64_t > get_first_derivative(const CLabels *lab, SGVector< float64_t > func, const TParameter *param) const
index_t vlen
Definition: SGVector.h:494
virtual SGVector< float64_t > get_predictive_variances(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab=NULL) const
#define ASSERT(x)
Definition: SGIO.h:201
virtual SGVector< float64_t > get_predictive_means(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab=NULL) const
double float64_t
Definition: common.h:50
virtual float64_t get_first_moment(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab, index_t i) const
static CGaussianLikelihood * obtain_from_generic(CLikelihoodModel *lik)
virtual SGVector< float64_t > get_second_derivative(const CLabels *lab, SGVector< float64_t > func, const TParameter *param) const
virtual SGVector< float64_t > get_log_zeroth_moments(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab) const
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
#define SG_SERROR(...)
Definition: SGIO.h:179
static float64_t exp(float64_t x)
Definition: Math.h:621
virtual SGVector< float64_t > get_log_probability_derivative_f(const CLabels *lab, SGVector< float64_t > func, index_t i) const
#define SG_ADD(...)
Definition: SGObject.h:81
virtual SGVector< float64_t > get_third_derivative(const CLabels *lab, SGVector< float64_t > func, const TParameter *param) const
The Likelihood model base class.
void set_const(T const_elem)
Definition: SGVector.cpp:152
static const float64_t PI
Definition: Math.h:2055

SHOGUN Machine Learning Toolbox - Documentation