49 using namespace Eigen;
54 CLogitVGPiecewiseBoundLikelihood::CLogitVGPiecewiseBoundLikelihood()
75 m_bound(0,0)=0.000188712193000;
76 m_bound(1,0)=0.028090310300000;
77 m_bound(2,0)=0.110211757000000;
78 m_bound(3,0)=0.232736440000000;
79 m_bound(4,0)=0.372524706000000;
80 m_bound(5,0)=0.504567936000000;
81 m_bound(6,0)=0.606280283000000;
82 m_bound(7,0)=0.666125432000000;
83 m_bound(8,0)=0.689334264000000;
84 m_bound(9,0)=0.693147181000000;
85 m_bound(10,0)=0.693147181000000;
86 m_bound(11,0)=0.689334264000000;
87 m_bound(12,0)=0.666125432000000;
88 m_bound(13,0)=0.606280283000000;
89 m_bound(14,0)=0.504567936000000;
90 m_bound(15,0)=0.372524706000000;
91 m_bound(16,0)=0.232736440000000;
92 m_bound(17,0)=0.110211757000000;
93 m_bound(18,0)=0.028090310400000;
94 m_bound(19,0)=0.000188712000000;
97 m_bound(1,1)=0.006648614600000;
98 m_bound(2,1)=0.034432684600000;
99 m_bound(3,1)=0.088701969900000;
100 m_bound(4,1)=0.168024214000000;
101 m_bound(5,1)=0.264032863000000;
102 m_bound(6,1)=0.360755794000000;
103 m_bound(7,1)=0.439094482000000;
104 m_bound(8,1)=0.485091758000000;
105 m_bound(9,1)=0.499419205000000;
106 m_bound(10,1)=0.500580795000000;
107 m_bound(11,1)=0.514908242000000;
108 m_bound(12,1)=0.560905518000000;
109 m_bound(13,1)=0.639244206000000;
110 m_bound(14,1)=0.735967137000000;
111 m_bound(15,1)=0.831975786000000;
112 m_bound(16,1)=0.911298030000000;
113 m_bound(17,1)=0.965567315000000;
114 m_bound(18,1)=0.993351385000000;
115 m_bound(19,1)=1.000000000000000;
118 m_bound(1,2)=0.000397791059000;
119 m_bound(2,2)=0.002753100850000;
120 m_bound(3,2)=0.008770186980000;
121 m_bound(4,2)=0.020034759300000;
122 m_bound(5,2)=0.037511596000000;
123 m_bound(6,2)=0.060543032900000;
124 m_bound(7,2)=0.086256780600000;
125 m_bound(8,2)=0.109213531000000;
126 m_bound(9,2)=0.123026104000000;
127 m_bound(10,2)=0.123026104000000;
128 m_bound(11,2)=0.109213531000000;
129 m_bound(12,2)=0.086256780600000;
130 m_bound(13,2)=0.060543032900000;
131 m_bound(14,2)=0.037511596000000;
132 m_bound(15,2)=0.020034759300000;
133 m_bound(16,2)=0.008770186980000;
134 m_bound(17,2)=0.002753100850000;
135 m_bound(18,2)=0.000397791059000;
139 m_bound(1,3)=-8.575194939999999;
140 m_bound(2,3)=-5.933689180000000;
141 m_bound(3,3)=-4.525933600000000;
142 m_bound(4,3)=-3.528107790000000;
143 m_bound(5,3)=-2.751548540000000;
144 m_bound(6,3)=-2.097898790000000;
145 m_bound(7,3)=-1.519690830000000;
146 m_bound(8,3)=-0.989533382000000;
147 m_bound(9,3)=-0.491473077000000;
149 m_bound(11,3)=0.491473077000000;
150 m_bound(12,3)=0.989533382000000;
151 m_bound(13,3)=1.519690830000000;
152 m_bound(14,3)=2.097898790000000;
153 m_bound(15,3)=2.751548540000000;
154 m_bound(16,3)=3.528107790000000;
155 m_bound(17,3)=4.525933600000000;
156 m_bound(18,3)=5.933689180000000;
157 m_bound(19,3)=8.575194939999999;
159 m_bound(0,4)=-8.575194939999999;
160 m_bound(1,4)=-5.933689180000000;
161 m_bound(2,4)=-4.525933600000000;
162 m_bound(3,4)=-3.528107790000000;
163 m_bound(4,4)=-2.751548540000000;
164 m_bound(5,4)=-2.097898790000000;
165 m_bound(6,4)=-1.519690830000000;
166 m_bound(7,4)=-0.989533382000000;
167 m_bound(8,4)=-0.491473077000000;
169 m_bound(10,4)=0.491473077000000;
170 m_bound(11,4)=0.989533382000000;
171 m_bound(12,4)=1.519690830000000;
172 m_bound(13,4)=2.097898790000000;
173 m_bound(14,4)=2.751548540000000;
174 m_bound(15,4)=3.528107790000000;
175 m_bound(16,4)=4.525933600000000;
176 m_bound(17,4)=5.933689180000000;
177 m_bound(18,4)=8.575194939999999;
199 const Map<MatrixXd> eigen_cdf_diff(m_cdf_diff.
matrix, num_rows, num_cols);
200 const Map<MatrixXd> eigen_pl(m_pl.
matrix, num_rows, num_cols);
201 const Map<MatrixXd> eigen_ph(m_ph.
matrix, num_rows, num_cols);
207 float64_t h_bak = eigen_h(eigen_h.size()-1);
209 eigen_h(eigen_h.size()-1) = 0;
212 const Map<MatrixXd> & eigen_ex0 = eigen_cdf_diff;
215 MatrixXd eigen_ex1 = ((eigen_pl - eigen_ph).array().rowwise()*eigen_s2.array().transpose()
216 + eigen_cdf_diff.array().rowwise()*eigen_mu.array().transpose()).matrix();
221 MatrixXd eigen_ex2 = ((eigen_mu.replicate(1,eigen_l.rows()).array().transpose().colwise() + eigen_l.array())*eigen_pl.array()
222 - (eigen_mu.replicate(1,eigen_h.rows()).array().transpose().colwise() + eigen_h.array())*eigen_ph.array()).matrix();
225 eigen_ex2 = (eigen_ex2.array().rowwise()*eigen_s2.array().transpose()).matrix();
228 eigen_ex2 += (eigen_cdf_diff.array().rowwise()*(eigen_mu.array().pow(2)
229 + eigen_s2.array()).transpose()).matrix();
235 eigen_f = (eigen_ex2.array().colwise()*eigen_a.array()
236 + eigen_ex1.array().colwise()*eigen_b.array()
237 + eigen_ex0.array().colwise()*eigen_c.array()).colwise().sum().matrix();
240 eigen_h(eigen_h.size()-1) = h_bak;
243 eigen_f = eigen_lab.cwiseProduct(eigen_mu) - eigen_f;
254 REQUIRE(param,
"Param is required (param should not be NULL)\n");
255 REQUIRE(param->
m_name,
"Param name is required (param->m_name should not be NULL)\n");
258 "Can't compute derivative of the variational expection ",
259 "of log LogitLikelihood using the piecewise bound ",
260 "wrt %s.%s parameter. The function only accepts mu and sigma2 as parameter",
275 const Map<MatrixXd> eigen_cdf_diff(m_cdf_diff.
matrix, num_rows, num_cols);
276 const Map<MatrixXd> eigen_pl(m_pl.
matrix, num_rows, num_cols);
277 const Map<MatrixXd> eigen_ph(m_ph.
matrix, num_rows, num_cols);
278 const Map<MatrixXd> eigen_weighted_pdf_diff(m_weighted_pdf_diff.
matrix, num_rows, num_cols);
279 const Map<MatrixXd> eigen_h2_plus_s2(m_h2_plus_s2.
matrix, num_rows, num_cols);
280 const Map<MatrixXd> eigen_l2_plus_s2(m_l2_plus_s2.
matrix, num_rows, num_cols);
285 float64_t h_bak = eigen_h(eigen_h.size()-1);
287 eigen_h(eigen_h.size()-1) = 0;
291 if (strcmp(param->
m_name,
"mu") == 0)
296 MatrixXd eigen_dmu2 = ((eigen_l2_plus_s2.array().rowwise()+eigen_s2.array().transpose())*eigen_pl.array()
297 - (eigen_h2_plus_s2.array().rowwise()+eigen_s2.array().transpose())*eigen_ph.array()).matrix();
299 eigen_dmu2 += (eigen_cdf_diff.array().rowwise()*(2.0*eigen_mu).array().transpose()).matrix();
302 Map<VectorXd> eigen_gmu(gmu.
vector, gmu.
vlen);
306 eigen_gmu = ((eigen_dmu2.array().colwise()*eigen_a.array())
307 + ((eigen_weighted_pdf_diff + eigen_cdf_diff).array().colwise()*eigen_b.array())
308 + ( (eigen_pl - eigen_ph).array().colwise()*eigen_c.array())).colwise().sum().matrix();
311 eigen_gmu = eigen_lab - eigen_gmu;
318 MatrixXd eigen_gs2_0 = (((-eigen_mu).replicate(1,eigen_l.rows()).array().transpose().colwise() + eigen_l.array())*eigen_pl.array()
319 - ((-eigen_mu).replicate(1,eigen_h.rows()).array().transpose().colwise() + eigen_h.array())*eigen_ph.array()).matrix();
321 eigen_gs2_0 = (eigen_gs2_0.array().colwise()*eigen_c.array()).matrix();
324 MatrixXd tmpl = (eigen_l2_plus_s2 - (eigen_mu.replicate(1,eigen_l.rows()).array().transpose().colwise()*eigen_l.array()).matrix()
325 ).cwiseProduct(eigen_pl);
328 MatrixXd tmph = (eigen_h2_plus_s2 - (eigen_mu.replicate(1,eigen_h.rows()).array().transpose().colwise()*eigen_h.array()).matrix()
329 ).cwiseProduct(eigen_ph);
332 MatrixXd eigen_gs2_1 = ((tmpl - tmph).array().colwise()*eigen_b.array()).matrix();
335 MatrixXd eigen_gs2_2 = (tmpl.array().colwise()*eigen_l.array() - tmph.array().colwise()*eigen_h.array()).matrix();
338 eigen_gs2_2 = (eigen_gs2_2.array().colwise()*eigen_a.array()).matrix();
341 Map<VectorXd> eigen_gs2(gs2.
vector, gs2.
vlen);
345 eigen_gs2 = ((eigen_cdf_diff + 0.5*eigen_weighted_pdf_diff).array().colwise()*eigen_a.array()
346 + (eigen_gs2_0 + eigen_gs2_1 + eigen_gs2_2).array().rowwise()/(2.0*eigen_s2).array().transpose()
347 ).colwise().sum().matrix();
350 eigen_gs2 = -eigen_gs2;
353 eigen_h(eigen_h.size()-1) = h_bak;
364 "Labels must be type of CBinaryLabels\n");
371 for(
index_t i = 0; i < m_lab.size(); ++i)
382 void CLogitVGPiecewiseBoundLikelihood::init()
385 "Variational piecewise bound for logit likelihood",
388 "The pdf given the lower range and parameters(mu and variance)",
391 "The pdf given the higher range and parameters(mu and variance)",
393 SG_ADD(&m_cdf_diff,
"cdf_h_minus_cdf_l",
394 "The CDF difference between the lower and higher range given the parameters(mu and variance)",
396 SG_ADD(&m_l2_plus_s2,
"l2_plus_sigma2",
397 "The result of l^2 + sigma^2",
399 SG_ADD(&m_h2_plus_s2,
"h2_plus_sigma2",
400 "The result of h^2 + sigma^2",
402 SG_ADD(&m_weighted_pdf_diff,
"weighted_pdf_diff",
403 "The result of l*pdf(l_norm)-h*pdf(h_norm)",
409 void CLogitVGPiecewiseBoundLikelihood::precompute()
427 m_pl = SGMatrix<float64_t>(num_rows,num_cols);
428 m_ph = SGMatrix<float64_t>(num_rows,num_cols);
429 m_cdf_diff = SGMatrix<float64_t>(num_rows,num_cols);
430 m_l2_plus_s2 = SGMatrix<float64_t>(num_rows,num_cols);
431 m_h2_plus_s2 = SGMatrix<float64_t>(num_rows,num_cols);
432 m_weighted_pdf_diff = SGMatrix<float64_t>(num_rows,num_cols);
434 Map<MatrixXd> eigen_pl(m_pl.
matrix, num_rows, num_cols);
435 Map<MatrixXd> eigen_ph(m_ph.matrix, num_rows, num_cols);
436 Map<MatrixXd> eigen_cdf_diff(m_cdf_diff.matrix, num_rows, num_cols);
437 Map<MatrixXd> eigen_l2_plus_s2(m_l2_plus_s2.matrix, num_rows, num_cols);
438 Map<MatrixXd> eigen_h2_plus_s2(m_h2_plus_s2.matrix, num_rows, num_cols);
439 Map<MatrixXd> eigen_weighted_pdf_diff(m_weighted_pdf_diff.matrix, num_rows, num_cols);
441 SGMatrix<float64_t> zl(num_rows, num_cols);
442 Map<MatrixXd> eigen_zl(zl.matrix, num_rows, num_cols);
443 SGMatrix<float64_t> zh(num_rows, num_cols);
444 Map<MatrixXd> eigen_zh(zh.matrix, num_rows, num_cols);
447 eigen_zl = ((-eigen_mu).replicate(1,eigen_l.rows()).array().transpose().colwise() + eigen_l.array()).matrix();
449 eigen_zh = ((-eigen_mu).replicate(1,eigen_h.rows()).array().transpose().colwise() + eigen_h.array()).matrix();
451 VectorXd eigen_s_inv = eigen_s2.array().sqrt().inverse().matrix();
454 eigen_zl = (eigen_zl.array().rowwise()*eigen_s_inv.array().transpose()).matrix();
456 eigen_zh = (eigen_zh.array().rowwise()*eigen_s_inv.array().transpose()).matrix();
460 for (
index_t r = 0; r < zl.num_rows; r++)
462 for (
index_t c = 0; c < zl.num_cols; c++)
477 eigen_pl = (eigen_pl.array().rowwise()*eigen_s_inv.array().transpose()).matrix();
479 eigen_ph = (eigen_ph.array().rowwise()*eigen_s_inv.array().transpose()).matrix();
481 SGMatrix<float64_t> & cl = zl;
482 SGMatrix<float64_t> & ch = zh;
484 for (
index_t r = 0; r < zl.num_rows; r++)
486 for (
index_t c = 0; c < zl.num_cols; c++)
495 Map<MatrixXd> eigen_cl(cl.matrix, num_rows, num_cols);
496 Map<MatrixXd> eigen_ch(ch.matrix, num_rows, num_cols);
498 eigen_cdf_diff = eigen_ch - eigen_cl;
503 float64_t h_bak = eigen_h(eigen_h.size()-1);
504 eigen_h(eigen_h.size()-1) = 0;
507 eigen_l2_plus_s2 = (eigen_s2.replicate(1,eigen_l.rows()).array().transpose().colwise() + eigen_l.array().pow(2)).matrix();
509 eigen_h2_plus_s2 = (eigen_s2.replicate(1,eigen_h.rows()).array().transpose().colwise() + eigen_h.array().pow(2)).matrix();
511 eigen_weighted_pdf_diff = (eigen_pl.array().colwise() * eigen_l.array() - eigen_ph.array().colwise() * eigen_h.array()).matrix();
514 eigen_h(eigen_h.size()-1) = h_bak;
Class that models Logit likelihood.
virtual ELabelType get_label_type() const =0
The class Labels models labels, i.e. class assignments of objects.
static const float64_t INFTY
infinity
virtual CSGObject * clone()
The variational Gaussian Likelihood base class. The variational distribution is Gaussian.
void set_default_variational_bound()
virtual void set_variational_distribution(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab)
virtual void set_variational_bound(SGMatrix< float64_t > bound)
virtual void set_likelihood(CLikelihoodModel *lik)
virtual SGVector< float64_t > get_variational_first_derivative(const TParameter *param) const
SGVector< float64_t > m_lab
virtual SGVector< float64_t > get_first_derivative_wrt_hyperparameter(const TParameter *param) const
virtual SGVector< float64_t > get_variational_expection()
static T max(T a, T b)
return the maximum of two integers
static float64_t univariate_log_pdf(float64_t sample, float64_t mu=0.0, float64_t sigma2=1.0)
virtual void set_variational_distribution(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab)
virtual const char * get_name() const
static float64_t exp(float64_t x)
static float64_t normal_cdf(float64_t x, float64_t std_dev=1)
Binary Labels for binary classification.
virtual ~CLogitVGPiecewiseBoundLikelihood()
SGVector< float64_t > m_mu
virtual void init_likelihood()
SGVector< float64_t > m_s2
T * get_column_vector(index_t col) const
static T abs(T a)
return the absolute value of a number