SHOGUN
v2.0.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Copyright (C) 2012 Jacob Walker 00008 * 00009 * Adapted from the GPML toolbox, specifically likT.m 00010 * http://www.gaussianprocess.org/gpml/code/matlab/doc/ 00011 * 00012 */ 00013 00014 00015 #include <shogun/regression/gp/StudentsTLikelihood.h> 00016 #ifdef HAVE_EIGEN3 00017 #include <shogun/modelselection/ParameterCombination.h> 00018 #include <shogun/mathematics/Statistics.h> 00019 #include <shogun/base/Parameter.h> 00020 #include <shogun/mathematics/eigen3.h> 00021 00022 using namespace shogun; 00023 using namespace Eigen; 00024 00025 CStudentsTLikelihood::CStudentsTLikelihood() : CLikelihoodModel() 00026 { 00027 init(); 00028 } 00029 00030 void CStudentsTLikelihood::init() 00031 { 00032 m_df = 3.0; 00033 m_sigma = 0.01; 00034 SG_ADD(&m_sigma, "sigma", "Observation Noise.", MS_AVAILABLE); 00035 } 00036 00037 CStudentsTLikelihood::~CStudentsTLikelihood() 00038 { 00039 } 00040 00041 00042 SGVector<float64_t> CStudentsTLikelihood::evaluate_means( 00043 SGVector<float64_t>& means) 00044 { 00045 return SGVector<float64_t>(means); 00046 } 00047 00048 SGVector<float64_t> CStudentsTLikelihood::evaluate_variances( 00049 SGVector<float64_t>& vars) 00050 { 00051 SGVector<float64_t> result(vars); 00052 00053 for (index_t i = 0; i < result.vlen; i++) 00054 { 00055 if (m_df < 2) 00056 result[i] = CMath::INFTY; 00057 else 00058 result[i] += m_df*(m_sigma*m_sigma)/float64_t(m_df-2); 00059 00060 } 00061 00062 return result; 00063 } 00064 00065 float64_t CStudentsTLikelihood::get_log_probability_f( 00066 CRegressionLabels* labels, SGVector<float64_t> m_function) 00067 { 00068 Map<VectorXd> function(m_function.vector, m_function.vlen); 00069 00070 float64_t temp = lgamma(m_df/2.0+0.5) - 00071 lgamma(m_df/2.0) - log(m_df*CMath::PI*m_sigma*m_sigma)/2.0; 00072 00073 VectorXd result(function.rows()); 00074 00075 for (index_t i = 0; i < function.rows(); i++) 00076 result[i] = labels->get_labels()[i] - function[i]; 00077 00078 result = result.cwiseProduct(result); 00079 00080 result /= m_df*m_sigma*m_sigma; 00081 00082 for (index_t i = 0; i < function.rows(); i++) 00083 { 00084 result[i] = -(m_df+1)*log(1.0+float64_t(result[i]))/2.0; 00085 00086 result[i] += temp; 00087 } 00088 00089 return result.sum(); 00090 } 00091 00092 SGVector<float64_t> CStudentsTLikelihood::get_log_probability_derivative_f( 00093 CRegressionLabels* labels, SGVector<float64_t> m_function, index_t j) 00094 { 00095 Map<VectorXd> function(m_function.vector, m_function.vlen); 00096 VectorXd result(function.rows()); 00097 00098 for (index_t i = 0; i < function.rows(); i++) 00099 result[i] = (labels->get_labels()[i] - function[i]); 00100 00101 VectorXd result_squared = result.cwiseProduct(result); 00102 00103 VectorXd a(function.rows()); 00104 VectorXd b(function.rows()); 00105 VectorXd c(function.rows()); 00106 VectorXd d(function.rows()); 00107 00108 SGVector<float64_t> sgresult(result.rows()); 00109 00110 for (index_t i = 0; i < function.rows(); i++) 00111 a[i] = result_squared[i] + m_df*m_sigma*m_sigma; 00112 00113 if (j == 1) 00114 { 00115 result = (m_df+1)*result.cwiseQuotient(a); 00116 for (index_t i = 0; i < result.rows(); i++) 00117 sgresult[i] = result[i]; 00118 return sgresult; 00119 } 00120 00121 for (index_t i = 0; i < function.rows(); i++) 00122 b[i] = result_squared[i] - m_df*m_sigma*m_sigma; 00123 00124 if (j == 2) 00125 { 00126 result = (m_df+1)*b.cwiseQuotient(a.cwiseProduct(a)); 00127 for (index_t i = 0; i < result.rows(); i++) 00128 sgresult[i] = result[i]; 00129 return sgresult; 00130 } 00131 00132 for (index_t i = 0; i < function.rows(); i++) 00133 c[i] = result_squared[i] - 3*m_df*m_sigma*m_sigma; 00134 00135 d = a.cwiseProduct(a); 00136 00137 if (j == 3) 00138 { 00139 result = (m_df+1)*2*result.cwiseProduct(c).cwiseQuotient(d.cwiseProduct(a)); 00140 for (index_t i = 0; i < result.rows(); i++) 00141 sgresult[i] = result[i]; 00142 return sgresult; 00143 00144 } 00145 00146 else 00147 { 00148 SG_ERROR("Invalid index for derivative\n"); 00149 return sgresult; 00150 } 00151 } 00152 00153 //Taken in log space then converted back to direct derivative 00154 SGVector<float64_t> CStudentsTLikelihood::get_first_derivative( 00155 CRegressionLabels* labels, TParameter* param, 00156 CSGObject* obj, SGVector<float64_t> m_function) 00157 { 00158 Map<VectorXd> function(m_function.vector, m_function.vlen); 00159 00160 SGVector<float64_t> sgresult(function.rows()); 00161 00162 VectorXd result(function.rows()); 00163 00164 for (index_t i = 0; i < function.rows(); i++) 00165 result[i] = (labels->get_labels()[i] - function[i]); 00166 00167 00168 VectorXd result_squared = result.cwiseProduct(result); 00169 00170 VectorXd a(function.rows()); 00171 VectorXd b(function.rows()); 00172 VectorXd c(function.rows()); 00173 VectorXd d(function.rows()); 00174 00175 00176 if (strcmp(param->m_name, "df") == 0 && obj == this) 00177 { 00178 for (index_t i = 0; i < function.rows(); i++) 00179 a[i] = result_squared[i] + m_df*m_sigma*m_sigma; 00180 00181 a = result_squared.cwiseQuotient(a); 00182 00183 a *= (m_df/2.0+.5); 00184 00185 for (index_t i = 0; i < function.rows(); i++) 00186 a[i] += m_df*( CStatistics::dlgamma(m_df/2.0+1/2.0)- 00187 CStatistics::dlgamma(m_df/2.0) )/2.0 - 1/2.0 00188 -m_df*log(1+result_squared[i]/(m_df*m_sigma*m_sigma))/2.0; 00189 00190 a *= (1-1/m_df); 00191 00192 result = a/(m_df-1); 00193 00194 for (index_t i = 0; i < result.rows(); i++) 00195 sgresult[i] = result[i]; 00196 return sgresult; 00197 } 00198 00199 00200 if (strcmp(param->m_name, "sigma") == 0 && obj == this) 00201 { 00202 for (index_t i = 0; i < function.rows(); i++) 00203 a[i] = result_squared[i] + m_df*m_sigma*m_sigma; 00204 00205 a = (m_df+1)*result_squared.cwiseQuotient(a); 00206 00207 for (index_t i = 0; i < function.rows(); i++) 00208 a[i] -= 1.0; 00209 00210 result = a/(m_sigma); 00211 00212 for (index_t i = 0; i < result.rows(); i++) 00213 sgresult[i] = result[i]; 00214 00215 return sgresult; 00216 } 00217 00218 00219 sgresult[0] = CMath::INFTY; 00220 00221 return sgresult; 00222 } 00223 00224 //Taken in log space then converted back to direct derivative 00225 SGVector<float64_t> CStudentsTLikelihood::get_second_derivative( 00226 CRegressionLabels* labels, TParameter* param, 00227 CSGObject* obj, SGVector<float64_t> m_function) 00228 { 00229 Map<VectorXd> function(m_function.vector, m_function.vlen); 00230 00231 SGVector<float64_t> sgresult(function.rows()); 00232 VectorXd result(function.rows()); 00233 00234 for (index_t i = 0; i < function.rows(); i++) 00235 result[i] = (labels->get_labels()[i] - function[i]); 00236 00237 VectorXd result_squared = result.cwiseProduct(result); 00238 00239 VectorXd a(function.rows()); 00240 VectorXd b(function.rows()); 00241 VectorXd c(function.rows()); 00242 VectorXd d(function.rows()); 00243 00244 if (strcmp(param->m_name, "df") == 0 && obj == this) 00245 { 00246 for (index_t i = 0; i < function.rows(); i++) 00247 a[i] = result_squared[i] + m_df*m_sigma*m_sigma; 00248 00249 b = result_squared.cwiseProduct(result_squared); 00250 00251 b = b - 3*result_squared*m_sigma*m_sigma*(m_df+1); 00252 00253 for (index_t i = 0; i < function.rows(); i++) 00254 b[i] = result_squared[i] - 3*m_sigma*m_sigma*(1+m_df); 00255 00256 b = b.cwiseProduct(result_squared); 00257 00258 for (index_t i = 0; i < function.rows(); i++) 00259 b[i] = b[i] + pow(m_sigma, 4)*m_df; 00260 00261 b *= m_df; 00262 00263 c = a.cwiseProduct(a); 00264 00265 c = c.cwiseProduct(a); 00266 00267 result = b.cwiseQuotient(c); 00268 00269 result = result/(m_df-1); 00270 00271 for (index_t i = 0; i < result.rows(); i++) 00272 sgresult[i] = result[i]; 00273 return sgresult; 00274 } 00275 00276 if (strcmp(param->m_name, "sigma") == 0 && obj == this) 00277 { 00278 for (index_t i = 0; i < function.rows(); i++) 00279 a[i] = result_squared[i] + m_df*m_sigma*m_sigma; 00280 00281 c = a.cwiseProduct(a); 00282 00283 c = c.cwiseProduct(a); 00284 00285 for (index_t i = 0; i < function.rows(); i++) 00286 b[i] = m_df*m_sigma*m_sigma - 3*result_squared[i]; 00287 00288 b *= m_sigma*m_sigma*m_df*2.0*(m_df+1); 00289 00290 result = b.cwiseQuotient(c)/m_sigma; 00291 00292 for (index_t i = 0; i < result.rows(); i++) 00293 sgresult[i] = result[i]; 00294 00295 return sgresult; 00296 00297 } 00298 00299 00300 sgresult[0] = CMath::INFTY; 00301 return sgresult; 00302 00303 } 00304 00305 #endif //HAVE_EIGEN3 00306 00307 00308