SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ProbitLikelihood.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  * Written (W) 2013 Roman Votyakov
8  */
9 
11 
12 #ifdef HAVE_EIGEN3
13 
17 
18 using namespace shogun;
19 using namespace Eigen;
20 
22 {
23 }
24 
26 {
27 }
28 
30  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
31 {
32  SGVector<float64_t> lp=get_log_zeroth_moments(mu, s2, lab);
33  Map<VectorXd> eigen_lp(lp.vector, lp.vlen);
34 
36  Map<VectorXd> eigen_r(r.vector, r.vlen);
37 
38  // evaluate predictive mean: ymu=2*p-1
39  eigen_r=2.0*eigen_lp.array().exp()-1.0;
40 
41  return r;
42 }
43 
45  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
46 {
47  SGVector<float64_t> lp=get_log_zeroth_moments(mu, s2, lab);
48  Map<VectorXd> eigen_lp(lp.vector, lp.vlen);
49 
51  Map<VectorXd> eigen_r(r.vector, r.vlen);
52 
53  // evaluate predictive variance: ys2=1-(2*p-1).^2
54  eigen_r=1-(2.0*eigen_lp.array().exp()-1.0).square();
55 
56  return r;
57 }
58 
60  SGVector<float64_t> func) const
61 {
62  // check the parameters
63  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
65  "Labels must be type of CBinaryLabels\n")
66  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
67  "length of the function vector\n")
68 
69  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
70  Map<VectorXd> eigen_y(y.vector, y.vlen);
71 
72  Map<VectorXd> eigen_f(func.vector, func.vlen);
73 
74  SGVector<float64_t> r(func.vlen);
75  Map<VectorXd> eigen_r(r.vector, r.vlen);
76 
77  // compute log pobability: log(normal_cdf(f.*y))
78  eigen_r=eigen_y.cwiseProduct(eigen_f);
79 
80  for (index_t i=0; i<eigen_r.size(); i++)
81  eigen_r[i]=CStatistics::lnormal_cdf(eigen_r[i]);
82 
83  return r;
84 }
85 
87  const CLabels* lab, SGVector<float64_t> func, index_t i) const
88 {
89  // check the parameters
90  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
92  "Labels must be type of CBinaryLabels\n")
93  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
94  "length of the function vector\n")
95  REQUIRE(i>=1 && i<=3, "Index for derivative should be 1, 2 or 3\n")
96 
97  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
98  Map<VectorXd> eigen_y(y.vector, y.vlen);
99 
100  Map<VectorXd> eigen_f(func.vector, func.vlen);
101 
102  SGVector<float64_t> dlp(func.vlen);
103  Map<VectorXd> eigen_dlp(dlp.vector, dlp.vlen);
104 
105  VectorXd eigen_yf=eigen_y.cwiseProduct(eigen_f);
106 
107  for (index_t j=0; j<eigen_yf.size(); j++)
108  {
109  float64_t v = eigen_yf[j];
111  {
112  //dlp( id2) = abs(den./num) * sqrt(2/pi); % strictly positive first derivative
113  eigen_dlp[j]=CMath::sqrt(2.0/CMath::PI)
115  }
116  else
117  {
118  //dlp(~id2) = exp(-z(~id2).*z(~id2)/2-lp(~id2))/sqrt(2*pi); % safe computation
119  eigen_dlp[j]=CMath::exp(-v*v/2.0-CStatistics::lnormal_cdf(v))
120  /CMath::sqrt(2.0*CMath::PI);
121  }
122  }
123 
124  SGVector<float64_t> r(func.vlen);
125  Map<VectorXd> eigen_r(r.vector, r.vlen);
126 
127  // compute derivatives of log probability wrt f
128 
129  if (i==1)
130  eigen_r=eigen_dlp;
131  else
132  //d2lp = -dlp.*abs(z+dlp); % strictly negative second derivative
133  eigen_r=(-eigen_dlp.array()*((eigen_yf.array()+eigen_dlp.array()).abs().array())).matrix();
134 
135  if (i==3)
136  //d3lp = -d2lp.*abs(z+2*dlp)-dlp; % strictly positive third derivative
137  eigen_r=(-eigen_r.array()*((eigen_yf.array()+2.0*eigen_dlp.array()).abs().array())
138  -eigen_dlp.array()).matrix();
139 
140  if (i==1 || i==3)
141  {
142  //varargout{2} = y.*varargout{2}
143  //varargout{4} = y.*varargout{4}
144  eigen_r=(eigen_r.array()*eigen_y.array()).matrix();
145  }
146 
147  return r;
148 }
149 
151  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
152 {
154 
155  if (lab)
156  {
157  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
158  "Length of the vector of means (%d), length of the vector of "
159  "variances (%d) and number of labels (%d) should be the same\n",
160  mu.vlen, s2.vlen, lab->get_num_labels())
162  "Labels must be type of CBinaryLabels\n")
163 
164  y=((CBinaryLabels*)lab)->get_labels();
165  }
166  else
167  {
168  REQUIRE(mu.vlen==s2.vlen, "Length of the vector of means (%d) and "
169  "length of the vector of variances (%d) should be the same\n",
170  mu.vlen, s2.vlen)
171 
173  y.set_const(1.0);
174  }
175 
176  Map<VectorXd> eigen_y(y.vector, y.vlen);
177  Map<VectorXd> eigen_mu(mu.vector, mu.vlen);
178  Map<VectorXd> eigen_s2(s2.vector, s2.vlen);
179 
181  Map<VectorXd> eigen_r(r.vector, r.vlen);
182 
183  // compute: lp=log(normal_cdf((mu.*y)./sqrt(1+sigma^2)))
184  eigen_r=eigen_mu.array()*eigen_y.array()/((1.0+eigen_s2.array()).sqrt());
185 
186  for (index_t i=0; i<eigen_r.size(); i++)
187  eigen_r[i]=CStatistics::lnormal_cdf(eigen_r[i]);
188 
189  return r;
190 }
191 
193  SGVector<float64_t> s2, const CLabels *lab, index_t i) const
194 {
195  // check the parameters
196  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
197  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
198  "Length of the vector of means (%d), length of the vector of "
199  "variances (%d) and number of labels (%d) should be the same\n",
200  mu.vlen, s2.vlen, lab->get_num_labels())
201  REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
203  "Labels must be type of CBinaryLabels\n")
204 
205  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
206 
207  float64_t z=y[i]*mu[i]/CMath::sqrt(1.0+s2[i]);
208 
209  // compute ncdf=normal_cdf(z)
211 
212  // compute npdf=normal_pdf(z)=(1/sqrt(2*pi))*exp(-z.^2/2)
213  float64_t npdf=(1.0/CMath::sqrt(2.0*CMath::PI))*CMath::exp(-0.5*CMath::sq(z));
214 
215  // compute the 1st moment: E[x] = mu + (y*s2*N(z))/(Phi(z)*sqrt(1+s2))
216  float64_t Ex=mu[i]+(npdf/ncdf)*(y[i]*s2[i])/CMath::sqrt(1.0+s2[i]);
217 
218  return Ex;
219 }
220 
222  SGVector<float64_t> s2, const CLabels *lab, index_t i) const
223 {
224  // check the parameters
225  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
226  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
227  "Length of the vector of means (%d), length of the vector of "
228  "variances (%d) and number of labels (%d) should be the same\n",
229  mu.vlen, s2.vlen, lab->get_num_labels())
230  REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
232  "Labels must be type of CBinaryLabels\n")
233 
234  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
235 
236  float64_t z=y[i]*mu[i]/CMath::sqrt(1.0+s2[i]);
237 
238  // compute ncdf=normal_cdf(z)
240 
241  // compute npdf=normal_pdf(z)=(1/sqrt(2*pi))*exp(-z.^2/2)
242  float64_t npdf=(1.0/CMath::sqrt(2.0*CMath::PI))*CMath::exp(-0.5*CMath::sq(z));
243 
244  SGVector<float64_t> r(y.vlen);
245  Map<VectorXd> eigen_r(r.vector, r.vlen);
246 
247  // compute the 2nd moment:
248  // Var[x] = s2 - (s2^2*N(z))/((1+s2)*Phi(z))*(z+N(z)/Phi(z))
249  float64_t Var=s2[i]-(CMath::sq(s2[i])/(1.0+s2[i]))*(npdf/ncdf)*(z+(npdf/ncdf));
250 
251  return Var;
252 }
253 
254 #endif /* HAVE_EIGEN3 */

SHOGUN Machine Learning Toolbox - Documentation