SHOGUN  v2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
StudentsTLikelihood.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  * Copyright (C) 2012 Jacob Walker
8  *
9  * Adapted from the GPML toolbox, specifically likT.m
10  * http://www.gaussianprocess.org/gpml/code/matlab/doc/
11  *
12  */
13 
14 
16 #ifdef HAVE_EIGEN3
19 #include <shogun/base/Parameter.h>
21 
22 using namespace shogun;
23 using namespace Eigen;
24 
26 {
27  init();
28 }
29 
30 void CStudentsTLikelihood::init()
31 {
32  m_df = 3.0;
33  m_sigma = 0.01;
34  SG_ADD(&m_sigma, "sigma", "Observation Noise.", MS_AVAILABLE);
35 }
36 
38 {
39 }
40 
41 
43  SGVector<float64_t>& means)
44 {
45  return SGVector<float64_t>(means);
46 }
47 
49  SGVector<float64_t>& vars)
50 {
51  SGVector<float64_t> result(vars);
52 
53  for (index_t i = 0; i < result.vlen; i++)
54  {
55  if (m_df < 2)
56  result[i] = CMath::INFTY;
57  else
58  result[i] += m_df*(m_sigma*m_sigma)/float64_t(m_df-2);
59 
60  }
61 
62  return result;
63 }
64 
66  CRegressionLabels* labels, SGVector<float64_t> m_function)
67 {
68  Map<VectorXd> function(m_function.vector, m_function.vlen);
69 
70  float64_t temp = lgamma(m_df/2.0+0.5) -
71  lgamma(m_df/2.0) - log(m_df*CMath::PI*m_sigma*m_sigma)/2.0;
72 
73  VectorXd result(function.rows());
74 
75  for (index_t i = 0; i < function.rows(); i++)
76  result[i] = labels->get_labels()[i] - function[i];
77 
78  result = result.cwiseProduct(result);
79 
80  result /= m_df*m_sigma*m_sigma;
81 
82  for (index_t i = 0; i < function.rows(); i++)
83  {
84  result[i] = -(m_df+1)*log(1.0+float64_t(result[i]))/2.0;
85 
86  result[i] += temp;
87  }
88 
89  return result.sum();
90 }
91 
93  CRegressionLabels* labels, SGVector<float64_t> m_function, index_t j)
94 {
95  Map<VectorXd> function(m_function.vector, m_function.vlen);
96  VectorXd result(function.rows());
97 
98  for (index_t i = 0; i < function.rows(); i++)
99  result[i] = (labels->get_labels()[i] - function[i]);
100 
101  VectorXd result_squared = result.cwiseProduct(result);
102 
103  VectorXd a(function.rows());
104  VectorXd b(function.rows());
105  VectorXd c(function.rows());
106  VectorXd d(function.rows());
107 
108  SGVector<float64_t> sgresult(result.rows());
109 
110  for (index_t i = 0; i < function.rows(); i++)
111  a[i] = result_squared[i] + m_df*m_sigma*m_sigma;
112 
113  if (j == 1)
114  {
115  result = (m_df+1)*result.cwiseQuotient(a);
116  for (index_t i = 0; i < result.rows(); i++)
117  sgresult[i] = result[i];
118  return sgresult;
119  }
120 
121  for (index_t i = 0; i < function.rows(); i++)
122  b[i] = result_squared[i] - m_df*m_sigma*m_sigma;
123 
124  if (j == 2)
125  {
126  result = (m_df+1)*b.cwiseQuotient(a.cwiseProduct(a));
127  for (index_t i = 0; i < result.rows(); i++)
128  sgresult[i] = result[i];
129  return sgresult;
130  }
131 
132  for (index_t i = 0; i < function.rows(); i++)
133  c[i] = result_squared[i] - 3*m_df*m_sigma*m_sigma;
134 
135  d = a.cwiseProduct(a);
136 
137  if (j == 3)
138  {
139  result = (m_df+1)*2*result.cwiseProduct(c).cwiseQuotient(d.cwiseProduct(a));
140  for (index_t i = 0; i < result.rows(); i++)
141  sgresult[i] = result[i];
142  return sgresult;
143 
144  }
145 
146  else
147  {
148  SG_ERROR("Invalid index for derivative\n");
149  return sgresult;
150  }
151 }
152 
153 //Taken in log space then converted back to direct derivative
155  CRegressionLabels* labels, TParameter* param,
156  CSGObject* obj, SGVector<float64_t> m_function)
157 {
158  Map<VectorXd> function(m_function.vector, m_function.vlen);
159 
160  SGVector<float64_t> sgresult(function.rows());
161 
162  VectorXd result(function.rows());
163 
164  for (index_t i = 0; i < function.rows(); i++)
165  result[i] = (labels->get_labels()[i] - function[i]);
166 
167 
168  VectorXd result_squared = result.cwiseProduct(result);
169 
170  VectorXd a(function.rows());
171  VectorXd b(function.rows());
172  VectorXd c(function.rows());
173  VectorXd d(function.rows());
174 
175 
176  if (strcmp(param->m_name, "df") == 0 && obj == this)
177  {
178  for (index_t i = 0; i < function.rows(); i++)
179  a[i] = result_squared[i] + m_df*m_sigma*m_sigma;
180 
181  a = result_squared.cwiseQuotient(a);
182 
183  a *= (m_df/2.0+.5);
184 
185  for (index_t i = 0; i < function.rows(); i++)
186  a[i] += m_df*( CStatistics::dlgamma(m_df/2.0+1/2.0)-
187  CStatistics::dlgamma(m_df/2.0) )/2.0 - 1/2.0
188  -m_df*log(1+result_squared[i]/(m_df*m_sigma*m_sigma))/2.0;
189 
190  a *= (1-1/m_df);
191 
192  result = a/(m_df-1);
193 
194  for (index_t i = 0; i < result.rows(); i++)
195  sgresult[i] = result[i];
196  return sgresult;
197  }
198 
199 
200  if (strcmp(param->m_name, "sigma") == 0 && obj == this)
201  {
202  for (index_t i = 0; i < function.rows(); i++)
203  a[i] = result_squared[i] + m_df*m_sigma*m_sigma;
204 
205  a = (m_df+1)*result_squared.cwiseQuotient(a);
206 
207  for (index_t i = 0; i < function.rows(); i++)
208  a[i] -= 1.0;
209 
210  result = a/(m_sigma);
211 
212  for (index_t i = 0; i < result.rows(); i++)
213  sgresult[i] = result[i];
214 
215  return sgresult;
216  }
217 
218 
219  sgresult[0] = CMath::INFTY;
220 
221  return sgresult;
222 }
223 
224 //Taken in log space then converted back to direct derivative
226  CRegressionLabels* labels, TParameter* param,
227  CSGObject* obj, SGVector<float64_t> m_function)
228 {
229  Map<VectorXd> function(m_function.vector, m_function.vlen);
230 
231  SGVector<float64_t> sgresult(function.rows());
232  VectorXd result(function.rows());
233 
234  for (index_t i = 0; i < function.rows(); i++)
235  result[i] = (labels->get_labels()[i] - function[i]);
236 
237  VectorXd result_squared = result.cwiseProduct(result);
238 
239  VectorXd a(function.rows());
240  VectorXd b(function.rows());
241  VectorXd c(function.rows());
242  VectorXd d(function.rows());
243 
244  if (strcmp(param->m_name, "df") == 0 && obj == this)
245  {
246  for (index_t i = 0; i < function.rows(); i++)
247  a[i] = result_squared[i] + m_df*m_sigma*m_sigma;
248 
249  b = result_squared.cwiseProduct(result_squared);
250 
251  b = b - 3*result_squared*m_sigma*m_sigma*(m_df+1);
252 
253  for (index_t i = 0; i < function.rows(); i++)
254  b[i] = result_squared[i] - 3*m_sigma*m_sigma*(1+m_df);
255 
256  b = b.cwiseProduct(result_squared);
257 
258  for (index_t i = 0; i < function.rows(); i++)
259  b[i] = b[i] + pow(m_sigma, 4)*m_df;
260 
261  b *= m_df;
262 
263  c = a.cwiseProduct(a);
264 
265  c = c.cwiseProduct(a);
266 
267  result = b.cwiseQuotient(c);
268 
269  result = result/(m_df-1);
270 
271  for (index_t i = 0; i < result.rows(); i++)
272  sgresult[i] = result[i];
273  return sgresult;
274  }
275 
276  if (strcmp(param->m_name, "sigma") == 0 && obj == this)
277  {
278  for (index_t i = 0; i < function.rows(); i++)
279  a[i] = result_squared[i] + m_df*m_sigma*m_sigma;
280 
281  c = a.cwiseProduct(a);
282 
283  c = c.cwiseProduct(a);
284 
285  for (index_t i = 0; i < function.rows(); i++)
286  b[i] = m_df*m_sigma*m_sigma - 3*result_squared[i];
287 
288  b *= m_sigma*m_sigma*m_df*2.0*(m_df+1);
289 
290  result = b.cwiseQuotient(c)/m_sigma;
291 
292  for (index_t i = 0; i < result.rows(); i++)
293  sgresult[i] = result[i];
294 
295  return sgresult;
296 
297  }
298 
299 
300  sgresult[0] = CMath::INFTY;
301  return sgresult;
302 
303 }
304 
305 #endif //HAVE_EIGEN3
306 
307 
308 

SHOGUN Machine Learning Toolbox - Documentation