SHOGUN
v3.2.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 * Written (W) 2013 Roman Votyakov 00008 */ 00009 00010 #include <shogun/machine/gp/ProbitLikelihood.h> 00011 00012 #ifdef HAVE_EIGEN3 00013 00014 #include <shogun/labels/BinaryLabels.h> 00015 #include <shogun/mathematics/eigen3.h> 00016 #include <shogun/mathematics/Statistics.h> 00017 00018 using namespace shogun; 00019 using namespace Eigen; 00020 00021 CProbitLikelihood::CProbitLikelihood() 00022 { 00023 } 00024 00025 CProbitLikelihood::~CProbitLikelihood() 00026 { 00027 } 00028 00029 SGVector<float64_t> CProbitLikelihood::get_predictive_means( 00030 SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const 00031 { 00032 SGVector<float64_t> lp=get_log_zeroth_moments(mu, s2, lab); 00033 Map<VectorXd> eigen_lp(lp.vector, lp.vlen); 00034 00035 SGVector<float64_t> r(lp.vlen); 00036 Map<VectorXd> eigen_r(r.vector, r.vlen); 00037 00038 // evaluate predictive mean: ymu=2*p-1 00039 eigen_r=2.0*eigen_lp.array().exp()-1.0; 00040 00041 return r; 00042 } 00043 00044 SGVector<float64_t> CProbitLikelihood::get_predictive_variances( 00045 SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const 00046 { 00047 SGVector<float64_t> lp=get_log_zeroth_moments(mu, s2, lab); 00048 Map<VectorXd> eigen_lp(lp.vector, lp.vlen); 00049 00050 SGVector<float64_t> r(lp.vlen); 00051 Map<VectorXd> eigen_r(r.vector, r.vlen); 00052 00053 // evaluate predictive variance: ys2=1-(2*p-1).^2 00054 eigen_r=1-(2.0*eigen_lp.array().exp()-1.0).square(); 00055 00056 return r; 00057 } 00058 00059 SGVector<float64_t> CProbitLikelihood::get_log_probability_f(const CLabels* lab, 00060 SGVector<float64_t> func) const 00061 { 00062 // check the parameters 00063 REQUIRE(lab, "Labels are required (lab should not be NULL)\n") 00064 REQUIRE(lab->get_label_type()==LT_BINARY, 00065 "Labels must be type of CBinaryLabels\n") 00066 REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match " 00067 "length of the function vector\n") 00068 00069 SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels(); 00070 Map<VectorXd> eigen_y(y.vector, y.vlen); 00071 00072 Map<VectorXd> eigen_f(func.vector, func.vlen); 00073 00074 SGVector<float64_t> r(func.vlen); 00075 Map<VectorXd> eigen_r(r.vector, r.vlen); 00076 00077 // compute log pobability: log(normal_cdf(f.*y)) 00078 eigen_r=eigen_y.cwiseProduct(eigen_f); 00079 00080 for (index_t i=0; i<eigen_r.size(); i++) 00081 eigen_r[i]=CStatistics::lnormal_cdf(eigen_r[i]); 00082 00083 return r; 00084 } 00085 00086 SGVector<float64_t> CProbitLikelihood::get_log_probability_derivative_f( 00087 const CLabels* lab, SGVector<float64_t> func, index_t i) const 00088 { 00089 // check the parameters 00090 REQUIRE(lab, "Labels are required (lab should not be NULL)\n") 00091 REQUIRE(lab->get_label_type()==LT_BINARY, 00092 "Labels must be type of CBinaryLabels\n") 00093 REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match " 00094 "length of the function vector\n") 00095 REQUIRE(i>=1 && i<=3, "Index for derivative should be 1, 2 or 3\n") 00096 00097 SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels(); 00098 Map<VectorXd> eigen_y(y.vector, y.vlen); 00099 00100 Map<VectorXd> eigen_f(func.vector, func.vlen); 00101 00102 SGVector<float64_t> r(func.vlen); 00103 Map<VectorXd> eigen_r(r.vector, r.vlen); 00104 00105 // compute ncdf=normal_cdf(y.*f) 00106 VectorXd eigen_ncdf=eigen_y.cwiseProduct(eigen_f); 00107 00108 for (index_t j=0; j<eigen_ncdf.size(); j++) 00109 eigen_ncdf[j]=CStatistics::normal_cdf(eigen_ncdf[j]); 00110 00111 // compute npdf=normal_pdf(f)=(1/sqrt(2*pi))*exp(-f.^2/2) 00112 VectorXd eigen_npdf=(1.0/CMath::sqrt(2.0*CMath::PI))* 00113 (-0.5*eigen_f.array().square()).exp(); 00114 00115 // compute z=npdf/ncdf 00116 VectorXd eigen_z=eigen_npdf.cwiseQuotient(eigen_ncdf); 00117 00118 // compute derivatives of log probability wrt f 00119 if (i == 1) 00120 { 00121 // compute the first derivative: dlp=y*z 00122 eigen_r=eigen_y.cwiseProduct(eigen_z); 00123 } 00124 else if (i == 2) 00125 { 00126 // compute the second derivative: d2lp=-z.^2-y.*f.*z 00127 eigen_r=-eigen_z.array().square()-eigen_y.array()*eigen_f.array()* 00128 eigen_z.array(); 00129 } 00130 else if (i == 3) 00131 { 00132 VectorXd eigen_z2=eigen_z.cwiseProduct(eigen_z); 00133 VectorXd eigen_z3=eigen_z2.cwiseProduct(eigen_z); 00134 00135 // compute the third derivative: d3lp=2*y.*z.^3+3*f.*z.^2+z.*y.*(f.^2-1) 00136 eigen_r=2.0*eigen_y.array()*eigen_z3.array()+3.0*eigen_f.array()* 00137 eigen_z2.array()+eigen_z.array()*eigen_y.array()* 00138 (eigen_f.array().square()-1.0); 00139 } 00140 00141 return r; 00142 } 00143 00144 SGVector<float64_t> CProbitLikelihood::get_log_zeroth_moments( 00145 SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const 00146 { 00147 SGVector<float64_t> y; 00148 00149 if (lab) 00150 { 00151 REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()), 00152 "Length of the vector of means (%d), length of the vector of " 00153 "variances (%d) and number of labels (%d) should be the same\n", 00154 mu.vlen, s2.vlen, lab->get_num_labels()) 00155 REQUIRE(lab->get_label_type()==LT_BINARY, 00156 "Labels must be type of CBinaryLabels\n") 00157 00158 y=((CBinaryLabels*)lab)->get_labels(); 00159 } 00160 else 00161 { 00162 REQUIRE(mu.vlen==s2.vlen, "Length of the vector of means (%d) and " 00163 "length of the vector of variances (%d) should be the same\n", 00164 mu.vlen, s2.vlen) 00165 00166 y=SGVector<float64_t>(mu.vlen); 00167 y.set_const(1.0); 00168 } 00169 00170 Map<VectorXd> eigen_y(y.vector, y.vlen); 00171 Map<VectorXd> eigen_mu(mu.vector, mu.vlen); 00172 Map<VectorXd> eigen_s2(s2.vector, s2.vlen); 00173 00174 SGVector<float64_t> r(y.vlen); 00175 Map<VectorXd> eigen_r(r.vector, r.vlen); 00176 00177 // compute: lp=log(normal_cdf((mu.*y)./sqrt(1+sigma^2))) 00178 eigen_r=eigen_mu.array()*eigen_y.array()/((1.0+eigen_s2.array()).sqrt()); 00179 00180 for (index_t i=0; i<eigen_r.size(); i++) 00181 eigen_r[i]=CStatistics::lnormal_cdf(eigen_r[i]); 00182 00183 return r; 00184 } 00185 00186 float64_t CProbitLikelihood::get_first_moment(SGVector<float64_t> mu, 00187 SGVector<float64_t> s2, const CLabels *lab, index_t i) const 00188 { 00189 // check the parameters 00190 REQUIRE(lab, "Labels are required (lab should not be NULL)\n") 00191 REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()), 00192 "Length of the vector of means (%d), length of the vector of " 00193 "variances (%d) and number of labels (%d) should be the same\n", 00194 mu.vlen, s2.vlen, lab->get_num_labels()) 00195 REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i) 00196 REQUIRE(lab->get_label_type()==LT_BINARY, 00197 "Labels must be type of CBinaryLabels\n") 00198 00199 SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels(); 00200 00201 float64_t z=y[i]*mu[i]/CMath::sqrt(1.0+s2[i]); 00202 00203 // compute ncdf=normal_cdf(z) 00204 float64_t ncdf=CStatistics::normal_cdf(z); 00205 00206 // compute npdf=normal_pdf(z)=(1/sqrt(2*pi))*exp(-z.^2/2) 00207 float64_t npdf=(1.0/CMath::sqrt(2.0*CMath::PI))*CMath::exp(-0.5*CMath::sq(z)); 00208 00209 // compute the 1st moment: E[x] = mu + (y*s2*N(z))/(Phi(z)*sqrt(1+s2)) 00210 float64_t Ex=mu[i]+(npdf/ncdf)*(y[i]*s2[i])/CMath::sqrt(1.0+s2[i]); 00211 00212 return Ex; 00213 } 00214 00215 float64_t CProbitLikelihood::get_second_moment(SGVector<float64_t> mu, 00216 SGVector<float64_t> s2, const CLabels *lab, index_t i) const 00217 { 00218 // check the parameters 00219 REQUIRE(lab, "Labels are required (lab should not be NULL)\n") 00220 REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()), 00221 "Length of the vector of means (%d), length of the vector of " 00222 "variances (%d) and number of labels (%d) should be the same\n", 00223 mu.vlen, s2.vlen, lab->get_num_labels()) 00224 REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i) 00225 REQUIRE(lab->get_label_type()==LT_BINARY, 00226 "Labels must be type of CBinaryLabels\n") 00227 00228 SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels(); 00229 00230 float64_t z=y[i]*mu[i]/CMath::sqrt(1.0+s2[i]); 00231 00232 // compute ncdf=normal_cdf(z) 00233 float64_t ncdf=CStatistics::normal_cdf(z); 00234 00235 // compute npdf=normal_pdf(z)=(1/sqrt(2*pi))*exp(-z.^2/2) 00236 float64_t npdf=(1.0/CMath::sqrt(2.0*CMath::PI))*CMath::exp(-0.5*CMath::sq(z)); 00237 00238 SGVector<float64_t> r(y.vlen); 00239 Map<VectorXd> eigen_r(r.vector, r.vlen); 00240 00241 // compute the 2nd moment: 00242 // Var[x] = s2 - (s2^2*N(z))/((1+s2)*Phi(z))*(z+N(z)/Phi(z)) 00243 float64_t Var=s2[i]-(CMath::sq(s2[i])/(1.0+s2[i]))*(npdf/ncdf)*(z+(npdf/ncdf)); 00244 00245 return Var; 00246 } 00247 00248 #endif /* HAVE_EIGEN3 */