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> r(func.vlen);
103  Map<VectorXd> eigen_r(r.vector, r.vlen);
104 
105  // compute ncdf=normal_cdf(y.*f)
106  VectorXd eigen_ncdf=eigen_y.cwiseProduct(eigen_f);
107 
108  for (index_t j=0; j<eigen_ncdf.size(); j++)
109  eigen_ncdf[j]=CStatistics::normal_cdf(eigen_ncdf[j]);
110 
111  // compute npdf=normal_pdf(f)=(1/sqrt(2*pi))*exp(-f.^2/2)
112  VectorXd eigen_npdf=(1.0/CMath::sqrt(2.0*CMath::PI))*
113  (-0.5*eigen_f.array().square()).exp();
114 
115  // compute z=npdf/ncdf
116  VectorXd eigen_z=eigen_npdf.cwiseQuotient(eigen_ncdf);
117 
118  // compute derivatives of log probability wrt f
119  if (i == 1)
120  {
121  // compute the first derivative: dlp=y*z
122  eigen_r=eigen_y.cwiseProduct(eigen_z);
123  }
124  else if (i == 2)
125  {
126  // compute the second derivative: d2lp=-z.^2-y.*f.*z
127  eigen_r=-eigen_z.array().square()-eigen_y.array()*eigen_f.array()*
128  eigen_z.array();
129  }
130  else if (i == 3)
131  {
132  VectorXd eigen_z2=eigen_z.cwiseProduct(eigen_z);
133  VectorXd eigen_z3=eigen_z2.cwiseProduct(eigen_z);
134 
135  // compute the third derivative: d3lp=2*y.*z.^3+3*f.*z.^2+z.*y.*(f.^2-1)
136  eigen_r=2.0*eigen_y.array()*eigen_z3.array()+3.0*eigen_f.array()*
137  eigen_z2.array()+eigen_z.array()*eigen_y.array()*
138  (eigen_f.array().square()-1.0);
139  }
140 
141  return r;
142 }
143 
145  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
146 {
148 
149  if (lab)
150  {
151  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
152  "Length of the vector of means (%d), length of the vector of "
153  "variances (%d) and number of labels (%d) should be the same\n",
154  mu.vlen, s2.vlen, lab->get_num_labels())
156  "Labels must be type of CBinaryLabels\n")
157 
158  y=((CBinaryLabels*)lab)->get_labels();
159  }
160  else
161  {
162  REQUIRE(mu.vlen==s2.vlen, "Length of the vector of means (%d) and "
163  "length of the vector of variances (%d) should be the same\n",
164  mu.vlen, s2.vlen)
165 
167  y.set_const(1.0);
168  }
169 
170  Map<VectorXd> eigen_y(y.vector, y.vlen);
171  Map<VectorXd> eigen_mu(mu.vector, mu.vlen);
172  Map<VectorXd> eigen_s2(s2.vector, s2.vlen);
173 
175  Map<VectorXd> eigen_r(r.vector, r.vlen);
176 
177  // compute: lp=log(normal_cdf((mu.*y)./sqrt(1+sigma^2)))
178  eigen_r=eigen_mu.array()*eigen_y.array()/((1.0+eigen_s2.array()).sqrt());
179 
180  for (index_t i=0; i<eigen_r.size(); i++)
181  eigen_r[i]=CStatistics::lnormal_cdf(eigen_r[i]);
182 
183  return r;
184 }
185 
187  SGVector<float64_t> s2, const CLabels *lab, index_t i) const
188 {
189  // check the parameters
190  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
191  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
192  "Length of the vector of means (%d), length of the vector of "
193  "variances (%d) and number of labels (%d) should be the same\n",
194  mu.vlen, s2.vlen, lab->get_num_labels())
195  REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
197  "Labels must be type of CBinaryLabels\n")
198 
199  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
200 
201  float64_t z=y[i]*mu[i]/CMath::sqrt(1.0+s2[i]);
202 
203  // compute ncdf=normal_cdf(z)
205 
206  // compute npdf=normal_pdf(z)=(1/sqrt(2*pi))*exp(-z.^2/2)
207  float64_t npdf=(1.0/CMath::sqrt(2.0*CMath::PI))*CMath::exp(-0.5*CMath::sq(z));
208 
209  // compute the 1st moment: E[x] = mu + (y*s2*N(z))/(Phi(z)*sqrt(1+s2))
210  float64_t Ex=mu[i]+(npdf/ncdf)*(y[i]*s2[i])/CMath::sqrt(1.0+s2[i]);
211 
212  return Ex;
213 }
214 
216  SGVector<float64_t> s2, const CLabels *lab, index_t i) const
217 {
218  // check the parameters
219  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
220  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
221  "Length of the vector of means (%d), length of the vector of "
222  "variances (%d) and number of labels (%d) should be the same\n",
223  mu.vlen, s2.vlen, lab->get_num_labels())
224  REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
226  "Labels must be type of CBinaryLabels\n")
227 
228  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
229 
230  float64_t z=y[i]*mu[i]/CMath::sqrt(1.0+s2[i]);
231 
232  // compute ncdf=normal_cdf(z)
234 
235  // compute npdf=normal_pdf(z)=(1/sqrt(2*pi))*exp(-z.^2/2)
236  float64_t npdf=(1.0/CMath::sqrt(2.0*CMath::PI))*CMath::exp(-0.5*CMath::sq(z));
237 
238  SGVector<float64_t> r(y.vlen);
239  Map<VectorXd> eigen_r(r.vector, r.vlen);
240 
241  // compute the 2nd moment:
242  // Var[x] = s2 - (s2^2*N(z))/((1+s2)*Phi(z))*(z+N(z)/Phi(z))
243  float64_t Var=s2[i]-(CMath::sq(s2[i])/(1.0+s2[i]))*(npdf/ncdf)*(z+(npdf/ncdf));
244 
245  return Var;
246 }
247 
248 #endif /* HAVE_EIGEN3 */

SHOGUN Machine Learning Toolbox - Documentation