SHOGUN  4.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
ProbitLikelihood.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (c) The Shogun Machine Learning Toolbox
3  * Written (W) 2013 Roman Votyakov
4  * All rights reserved.
5  *
6  * Redistribution and use in source and binary forms, with or without
7  * modification, are permitted provided that the following conditions are met:
8  *
9  * 1. Redistributions of source code must retain the above copyright notice, this
10  * list of conditions and the following disclaimer.
11  * 2. Redistributions in binary form must reproduce the above copyright notice,
12  * this list of conditions and the following disclaimer in the documentation
13  * and/or other materials provided with the distribution.
14  *
15  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
19  * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25  *
26  * The views and conclusions contained in the software and documentation are those
27  * of the authors and should not be interpreted as representing official policies,
28  * either expressed or implied, of the Shogun Development Team.
29  *
30  */
32 
33 
37 
38 using namespace shogun;
39 using namespace Eigen;
40 
42 {
43 }
44 
46 {
47 }
48 
50  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
51 {
52  SGVector<float64_t> lp=get_log_zeroth_moments(mu, s2, lab);
53  Map<VectorXd> eigen_lp(lp.vector, lp.vlen);
54 
56  Map<VectorXd> eigen_r(r.vector, r.vlen);
57 
58  // evaluate predictive mean: ymu=2*p-1
59  eigen_r=2.0*eigen_lp.array().exp()-1.0;
60 
61  return r;
62 }
63 
65  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
66 {
67  SGVector<float64_t> lp=get_log_zeroth_moments(mu, s2, lab);
68  Map<VectorXd> eigen_lp(lp.vector, lp.vlen);
69 
71  Map<VectorXd> eigen_r(r.vector, r.vlen);
72 
73  // evaluate predictive variance: ys2=1-(2*p-1).^2
74  eigen_r=1-(2.0*eigen_lp.array().exp()-1.0).square();
75 
76  return r;
77 }
78 
80  SGVector<float64_t> func) const
81 {
82  // check the parameters
83  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
85  "Labels must be type of CBinaryLabels\n")
86  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
87  "length of the function vector\n")
88 
89  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
90  Map<VectorXd> eigen_y(y.vector, y.vlen);
91 
92  Map<VectorXd> eigen_f(func.vector, func.vlen);
93 
94  SGVector<float64_t> r(func.vlen);
95  Map<VectorXd> eigen_r(r.vector, r.vlen);
96 
97  // compute log pobability: log(normal_cdf(f.*y))
98  eigen_r=eigen_y.cwiseProduct(eigen_f);
99 
100  for (index_t i=0; i<eigen_r.size(); i++)
101  eigen_r[i]=CStatistics::lnormal_cdf(eigen_r[i]);
102 
103  return r;
104 }
105 
107  const CLabels* lab, SGVector<float64_t> func, index_t i) const
108 {
109  // check the parameters
110  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
112  "Labels must be type of CBinaryLabels\n")
113  REQUIRE(lab->get_num_labels()==func.vlen, "Number of labels must match "
114  "length of the function vector\n")
115  REQUIRE(i>=1 && i<=3, "Index for derivative should be 1, 2 or 3\n")
116 
117  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
118  Map<VectorXd> eigen_y(y.vector, y.vlen);
119 
120  Map<VectorXd> eigen_f(func.vector, func.vlen);
121 
122  SGVector<float64_t> dlp(func.vlen);
123  Map<VectorXd> eigen_dlp(dlp.vector, dlp.vlen);
124 
125  VectorXd eigen_yf=eigen_y.cwiseProduct(eigen_f);
126 
127  for (index_t j=0; j<eigen_yf.size(); j++)
128  {
129  float64_t v = eigen_yf[j];
131  {
132  //dlp( id2) = abs(den./num) * sqrt(2/pi); % strictly positive first derivative
133  eigen_dlp[j]=CMath::sqrt(2.0/CMath::PI)
135  }
136  else
137  {
138  //dlp(~id2) = exp(-z(~id2).*z(~id2)/2-lp(~id2))/sqrt(2*pi); % safe computation
139  eigen_dlp[j]=CMath::exp(-v*v/2.0-CStatistics::lnormal_cdf(v))
140  /CMath::sqrt(2.0*CMath::PI);
141  }
142  }
143 
144  SGVector<float64_t> r(func.vlen);
145  Map<VectorXd> eigen_r(r.vector, r.vlen);
146 
147  // compute derivatives of log probability wrt f
148 
149  if (i==1)
150  eigen_r=eigen_dlp;
151  else
152  //d2lp = -dlp.*abs(z+dlp); % strictly negative second derivative
153  eigen_r=(-eigen_dlp.array()*((eigen_yf.array()+eigen_dlp.array()).abs().array())).matrix();
154 
155  if (i==3)
156  //d3lp = -d2lp.*abs(z+2*dlp)-dlp; % strictly positive third derivative
157  eigen_r=(-eigen_r.array()*((eigen_yf.array()+2.0*eigen_dlp.array()).abs().array())
158  -eigen_dlp.array()).matrix();
159 
160  if (i==1 || i==3)
161  {
162  //varargout{2} = y.*varargout{2}
163  //varargout{4} = y.*varargout{4}
164  eigen_r=(eigen_r.array()*eigen_y.array()).matrix();
165  }
166 
167  return r;
168 }
169 
171  SGVector<float64_t> mu, SGVector<float64_t> s2, const CLabels* lab) const
172 {
174 
175  if (lab)
176  {
177  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
178  "Length of the vector of means (%d), length of the vector of "
179  "variances (%d) and number of labels (%d) should be the same\n",
180  mu.vlen, s2.vlen, lab->get_num_labels())
182  "Labels must be type of CBinaryLabels\n")
183 
184  y=((CBinaryLabels*)lab)->get_labels();
185  }
186  else
187  {
188  REQUIRE(mu.vlen==s2.vlen, "Length of the vector of means (%d) and "
189  "length of the vector of variances (%d) should be the same\n",
190  mu.vlen, s2.vlen)
191 
193  y.set_const(1.0);
194  }
195 
196  Map<VectorXd> eigen_y(y.vector, y.vlen);
197  Map<VectorXd> eigen_mu(mu.vector, mu.vlen);
198  Map<VectorXd> eigen_s2(s2.vector, s2.vlen);
199 
201  Map<VectorXd> eigen_r(r.vector, r.vlen);
202 
203  // compute: lp=log(normal_cdf((mu.*y)./sqrt(1+sigma^2)))
204  eigen_r=eigen_mu.array()*eigen_y.array()/((1.0+eigen_s2.array()).sqrt());
205 
206  for (index_t i=0; i<eigen_r.size(); i++)
207  eigen_r[i]=CStatistics::lnormal_cdf(eigen_r[i]);
208 
209  return r;
210 }
211 
213  SGVector<float64_t> s2, const CLabels *lab, index_t i) const
214 {
215  // check the parameters
216  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
217  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
218  "Length of the vector of means (%d), length of the vector of "
219  "variances (%d) and number of labels (%d) should be the same\n",
220  mu.vlen, s2.vlen, lab->get_num_labels())
221  REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
223  "Labels must be type of CBinaryLabels\n")
224 
225  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
226 
227  float64_t z=y[i]*mu[i]/CMath::sqrt(1.0+s2[i]);
228 
229  // compute ncdf=normal_cdf(z)
231 
232  // compute npdf=normal_pdf(z)=(1/sqrt(2*pi))*exp(-z.^2/2)
233  float64_t npdf=(1.0/CMath::sqrt(2.0*CMath::PI))*CMath::exp(-0.5*CMath::sq(z));
234 
235  // compute the 1st moment: E[x] = mu + (y*s2*N(z))/(Phi(z)*sqrt(1+s2))
236  float64_t Ex=mu[i]+(npdf/ncdf)*(y[i]*s2[i])/CMath::sqrt(1.0+s2[i]);
237 
238  return Ex;
239 }
240 
242  SGVector<float64_t> s2, const CLabels *lab, index_t i) const
243 {
244  // check the parameters
245  REQUIRE(lab, "Labels are required (lab should not be NULL)\n")
246  REQUIRE((mu.vlen==s2.vlen) && (mu.vlen==lab->get_num_labels()),
247  "Length of the vector of means (%d), length of the vector of "
248  "variances (%d) and number of labels (%d) should be the same\n",
249  mu.vlen, s2.vlen, lab->get_num_labels())
250  REQUIRE(i>=0 && i<=mu.vlen, "Index (%d) out of bounds!\n", i)
252  "Labels must be type of CBinaryLabels\n")
253 
254  SGVector<float64_t> y=((CBinaryLabels*)lab)->get_labels();
255 
256  float64_t z=y[i]*mu[i]/CMath::sqrt(1.0+s2[i]);
257 
258  // compute ncdf=normal_cdf(z)
260 
261  // compute npdf=normal_pdf(z)=(1/sqrt(2*pi))*exp(-z.^2/2)
262  float64_t npdf=(1.0/CMath::sqrt(2.0*CMath::PI))*CMath::exp(-0.5*CMath::sq(z));
263 
264  SGVector<float64_t> r(y.vlen);
265  Map<VectorXd> eigen_r(r.vector, r.vlen);
266 
267  // compute the 2nd moment:
268  // Var[x] = s2 - (s2^2*N(z))/((1+s2)*Phi(z))*(z+N(z)/Phi(z))
269  float64_t Var=s2[i]-(CMath::sq(s2[i])/(1.0+s2[i]))*(npdf/ncdf)*(z+(npdf/ncdf));
270 
271  return Var;
272 }
273 
virtual ELabelType get_label_type() const =0
binary labels +1/-1
Definition: LabelTypes.h:18
int32_t index_t
Definition: common.h:62
The class Labels models labels, i.e. class assignments of objects.
Definition: Labels.h:43
virtual int32_t get_num_labels() const =0
static T sq(T x)
Definition: Math.h:450
Definition: SGMatrix.h:20
#define REQUIRE(x,...)
Definition: SGIO.h:206
static const float64_t ERFC_CASE2
Definition: Statistics.h:408
virtual SGVector< float64_t > get_predictive_variances(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab=NULL) const
static float64_t lnormal_cdf(float64_t x)
Definition: Statistics.cpp:418
index_t vlen
Definition: SGVector.h:492
virtual SGVector< float64_t > get_log_probability_f(const CLabels *lab, SGVector< float64_t > func) const
double float64_t
Definition: common.h:50
virtual SGVector< float64_t > get_log_zeroth_moments(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab) const
virtual SGVector< float64_t > get_predictive_means(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab=NULL) const
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
static float64_t erfc8_weighted_sum(float64_t x)
Definition: Statistics.cpp:464
static float64_t exp(float64_t x)
Definition: Math.h:621
virtual float64_t get_first_moment(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab, index_t i) const
static float64_t normal_cdf(float64_t x, float64_t std_dev=1)
Definition: Statistics.cpp:508
Binary Labels for binary classification.
Definition: BinaryLabels.h:37
static float32_t sqrt(float32_t x)
Definition: Math.h:459
virtual SGVector< float64_t > get_log_probability_derivative_f(const CLabels *lab, SGVector< float64_t > func, index_t i) const
virtual float64_t get_second_moment(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab, index_t i) const
void set_const(T const_elem)
Definition: SGVector.cpp:150
static T abs(T a)
Definition: Math.h:179
static const float64_t PI
Definition: Math.h:2055

SHOGUN Machine Learning Toolbox - Documentation