50 using namespace Eigen;
55 CLogitVGPiecewiseBoundLikelihood::CLogitVGPiecewiseBoundLikelihood()
76 m_bound(0,0)=0.000188712193000;
77 m_bound(1,0)=0.028090310300000;
78 m_bound(2,0)=0.110211757000000;
79 m_bound(3,0)=0.232736440000000;
80 m_bound(4,0)=0.372524706000000;
81 m_bound(5,0)=0.504567936000000;
82 m_bound(6,0)=0.606280283000000;
83 m_bound(7,0)=0.666125432000000;
84 m_bound(8,0)=0.689334264000000;
85 m_bound(9,0)=0.693147181000000;
86 m_bound(10,0)=0.693147181000000;
87 m_bound(11,0)=0.689334264000000;
88 m_bound(12,0)=0.666125432000000;
89 m_bound(13,0)=0.606280283000000;
90 m_bound(14,0)=0.504567936000000;
91 m_bound(15,0)=0.372524706000000;
92 m_bound(16,0)=0.232736440000000;
93 m_bound(17,0)=0.110211757000000;
94 m_bound(18,0)=0.028090310400000;
95 m_bound(19,0)=0.000188712000000;
98 m_bound(1,1)=0.006648614600000;
99 m_bound(2,1)=0.034432684600000;
100 m_bound(3,1)=0.088701969900000;
101 m_bound(4,1)=0.168024214000000;
102 m_bound(5,1)=0.264032863000000;
103 m_bound(6,1)=0.360755794000000;
104 m_bound(7,1)=0.439094482000000;
105 m_bound(8,1)=0.485091758000000;
106 m_bound(9,1)=0.499419205000000;
107 m_bound(10,1)=0.500580795000000;
108 m_bound(11,1)=0.514908242000000;
109 m_bound(12,1)=0.560905518000000;
110 m_bound(13,1)=0.639244206000000;
111 m_bound(14,1)=0.735967137000000;
112 m_bound(15,1)=0.831975786000000;
113 m_bound(16,1)=0.911298030000000;
114 m_bound(17,1)=0.965567315000000;
115 m_bound(18,1)=0.993351385000000;
116 m_bound(19,1)=1.000000000000000;
119 m_bound(1,2)=0.000397791059000;
120 m_bound(2,2)=0.002753100850000;
121 m_bound(3,2)=0.008770186980000;
122 m_bound(4,2)=0.020034759300000;
123 m_bound(5,2)=0.037511596000000;
124 m_bound(6,2)=0.060543032900000;
125 m_bound(7,2)=0.086256780600000;
126 m_bound(8,2)=0.109213531000000;
127 m_bound(9,2)=0.123026104000000;
128 m_bound(10,2)=0.123026104000000;
129 m_bound(11,2)=0.109213531000000;
130 m_bound(12,2)=0.086256780600000;
131 m_bound(13,2)=0.060543032900000;
132 m_bound(14,2)=0.037511596000000;
133 m_bound(15,2)=0.020034759300000;
134 m_bound(16,2)=0.008770186980000;
135 m_bound(17,2)=0.002753100850000;
136 m_bound(18,2)=0.000397791059000;
140 m_bound(1,3)=-8.575194939999999;
141 m_bound(2,3)=-5.933689180000000;
142 m_bound(3,3)=-4.525933600000000;
143 m_bound(4,3)=-3.528107790000000;
144 m_bound(5,3)=-2.751548540000000;
145 m_bound(6,3)=-2.097898790000000;
146 m_bound(7,3)=-1.519690830000000;
147 m_bound(8,3)=-0.989533382000000;
148 m_bound(9,3)=-0.491473077000000;
150 m_bound(11,3)=0.491473077000000;
151 m_bound(12,3)=0.989533382000000;
152 m_bound(13,3)=1.519690830000000;
153 m_bound(14,3)=2.097898790000000;
154 m_bound(15,3)=2.751548540000000;
155 m_bound(16,3)=3.528107790000000;
156 m_bound(17,3)=4.525933600000000;
157 m_bound(18,3)=5.933689180000000;
158 m_bound(19,3)=8.575194939999999;
160 m_bound(0,4)=-8.575194939999999;
161 m_bound(1,4)=-5.933689180000000;
162 m_bound(2,4)=-4.525933600000000;
163 m_bound(3,4)=-3.528107790000000;
164 m_bound(4,4)=-2.751548540000000;
165 m_bound(5,4)=-2.097898790000000;
166 m_bound(6,4)=-1.519690830000000;
167 m_bound(7,4)=-0.989533382000000;
168 m_bound(8,4)=-0.491473077000000;
170 m_bound(10,4)=0.491473077000000;
171 m_bound(11,4)=0.989533382000000;
172 m_bound(12,4)=1.519690830000000;
173 m_bound(13,4)=2.097898790000000;
174 m_bound(14,4)=2.751548540000000;
175 m_bound(15,4)=3.528107790000000;
176 m_bound(16,4)=4.525933600000000;
177 m_bound(17,4)=5.933689180000000;
178 m_bound(18,4)=8.575194939999999;
208 float64_t h_bak = eigen_h(eigen_h.size()-1);
210 eigen_h(eigen_h.size()-1) = 0;
216 MatrixXd eigen_ex1 = ((eigen_pl - eigen_ph).array().rowwise()*eigen_s2.array().transpose()
217 + eigen_cdf_diff.array().rowwise()*eigen_mu.array().transpose()).matrix();
222 MatrixXd eigen_ex2 = ((eigen_mu.replicate(1,eigen_l.rows()).array().transpose().colwise() + eigen_l.array())*eigen_pl.array()
223 - (eigen_mu.replicate(1,eigen_h.rows()).array().transpose().colwise() + eigen_h.array())*eigen_ph.array()).matrix();
226 eigen_ex2 = (eigen_ex2.array().rowwise()*eigen_s2.array().transpose()).matrix();
229 eigen_ex2 += (eigen_cdf_diff.array().rowwise()*(eigen_mu.array().pow(2)
230 + eigen_s2.array()).transpose()).matrix();
236 eigen_f = (eigen_ex2.array().colwise()*eigen_a.array()
237 + eigen_ex1.array().colwise()*eigen_b.array()
238 + eigen_ex0.array().colwise()*eigen_c.array()).colwise().sum().matrix();
241 eigen_h(eigen_h.size()-1) = h_bak;
244 eigen_f = eigen_lab.cwiseProduct(eigen_mu) - eigen_f;
255 REQUIRE(param,
"Param is required (param should not be NULL)\n");
256 REQUIRE(param->
m_name,
"Param name is required (param->m_name should not be NULL)\n");
259 "Can't compute derivative of the variational expection ",
260 "of log LogitLikelihood using the piecewise bound ",
261 "wrt %s.%s parameter. The function only accepts mu and sigma2 as parameter",
279 const Map<MatrixXd> eigen_weighted_pdf_diff(m_weighted_pdf_diff.
matrix, num_rows, num_cols);
286 float64_t h_bak = eigen_h(eigen_h.size()-1);
288 eigen_h(eigen_h.size()-1) = 0;
292 if (strcmp(param->
m_name,
"mu") == 0)
297 MatrixXd eigen_dmu2 = ((eigen_l2_plus_s2.array().rowwise()+eigen_s2.array().transpose())*eigen_pl.array()
298 - (eigen_h2_plus_s2.array().rowwise()+eigen_s2.array().transpose())*eigen_ph.array()).matrix();
300 eigen_dmu2 += (eigen_cdf_diff.array().rowwise()*(2.0*eigen_mu).array().transpose()).matrix();
307 eigen_gmu = ((eigen_dmu2.array().colwise()*eigen_a.array())
308 + ((eigen_weighted_pdf_diff + eigen_cdf_diff).array().colwise()*eigen_b.array())
309 + ( (eigen_pl - eigen_ph).array().colwise()*eigen_c.array())).colwise().sum().matrix();
312 eigen_gmu = eigen_lab - eigen_gmu;
319 MatrixXd eigen_gs2_0 = (((-eigen_mu).replicate(1,eigen_l.rows()).array().transpose().colwise() + eigen_l.array())*eigen_pl.array()
320 - ((-eigen_mu).replicate(1,eigen_h.rows()).array().transpose().colwise() + eigen_h.array())*eigen_ph.array()).matrix();
322 eigen_gs2_0 = (eigen_gs2_0.array().colwise()*eigen_c.array()).matrix();
325 MatrixXd tmpl = (eigen_l2_plus_s2 - (eigen_mu.replicate(1,eigen_l.rows()).array().transpose().colwise()*eigen_l.array()).matrix()
326 ).cwiseProduct(eigen_pl);
329 MatrixXd tmph = (eigen_h2_plus_s2 - (eigen_mu.replicate(1,eigen_h.rows()).array().transpose().colwise()*eigen_h.array()).matrix()
330 ).cwiseProduct(eigen_ph);
333 MatrixXd eigen_gs2_1 = ((tmpl - tmph).array().colwise()*eigen_b.array()).matrix();
336 MatrixXd eigen_gs2_2 = (tmpl.array().colwise()*eigen_l.array() - tmph.array().colwise()*eigen_h.array()).matrix();
339 eigen_gs2_2 = (eigen_gs2_2.array().colwise()*eigen_a.array()).matrix();
346 eigen_gs2 = ((eigen_cdf_diff + 0.5*eigen_weighted_pdf_diff).array().colwise()*eigen_a.array()
347 + (eigen_gs2_0 + eigen_gs2_1 + eigen_gs2_2).array().rowwise()/(2.0*eigen_s2).array().transpose()
348 ).colwise().sum().matrix();
351 eigen_gs2 = -eigen_gs2;
354 eigen_h(eigen_h.size()-1) = h_bak;
367 "Labels must be type of CBinaryLabels\n");
374 for(
index_t i = 0; i < m_lab.size(); ++i)
388 void CLogitVGPiecewiseBoundLikelihood::init()
391 "Variational piecewise bound for logit likelihood",
394 "The pdf given the lower range and parameters(mu and variance)",
397 "The pdf given the higher range and parameters(mu and variance)",
399 SG_ADD(&m_cdf_diff,
"cdf_h_minus_cdf_l",
400 "The CDF difference between the lower and higher range given the parameters(mu and variance)",
402 SG_ADD(&m_l2_plus_s2,
"l2_plus_sigma2",
403 "The result of l^2 + sigma^2",
405 SG_ADD(&m_h2_plus_s2,
"h2_plus_sigma2",
406 "The result of h^2 + sigma^2",
408 SG_ADD(&m_weighted_pdf_diff,
"weighted_pdf_diff",
409 "The result of l*pdf(l_norm)-h*pdf(h_norm)",
415 void CLogitVGPiecewiseBoundLikelihood::precompute()
433 m_pl = SGMatrix<float64_t>(num_rows,num_cols);
434 m_ph = SGMatrix<float64_t>(num_rows,num_cols);
435 m_cdf_diff = SGMatrix<float64_t>(num_rows,num_cols);
436 m_l2_plus_s2 = SGMatrix<float64_t>(num_rows,num_cols);
437 m_h2_plus_s2 = SGMatrix<float64_t>(num_rows,num_cols);
438 m_weighted_pdf_diff = SGMatrix<float64_t>(num_rows,num_cols);
442 Map<MatrixXd> eigen_cdf_diff(m_cdf_diff.matrix, num_rows, num_cols);
443 Map<MatrixXd> eigen_l2_plus_s2(m_l2_plus_s2.matrix, num_rows, num_cols);
444 Map<MatrixXd> eigen_h2_plus_s2(m_h2_plus_s2.matrix, num_rows, num_cols);
445 Map<MatrixXd> eigen_weighted_pdf_diff(m_weighted_pdf_diff.matrix, num_rows, num_cols);
447 SGMatrix<float64_t> zl(num_rows, num_cols);
449 SGMatrix<float64_t> zh(num_rows, num_cols);
453 eigen_zl = ((-eigen_mu).replicate(1,eigen_l.rows()).array().transpose().colwise() + eigen_l.array()).matrix();
455 eigen_zh = ((-eigen_mu).replicate(1,eigen_h.rows()).array().transpose().colwise() + eigen_h.array()).matrix();
457 VectorXd eigen_s_inv = eigen_s2.array().sqrt().inverse().matrix();
460 eigen_zl = (eigen_zl.array().rowwise()*eigen_s_inv.array().transpose()).matrix();
462 eigen_zh = (eigen_zh.array().rowwise()*eigen_s_inv.array().transpose()).matrix();
466 for (
index_t r = 0; r < zl.num_rows; r++)
468 for (
index_t c = 0; c < zl.num_cols; c++)
483 eigen_pl = (eigen_pl.array().rowwise()*eigen_s_inv.array().transpose()).matrix();
485 eigen_ph = (eigen_ph.array().rowwise()*eigen_s_inv.array().transpose()).matrix();
487 SGMatrix<float64_t> & cl = zl;
488 SGMatrix<float64_t> & ch = zh;
490 for (
index_t r = 0; r < zl.num_rows; r++)
492 for (
index_t c = 0; c < zl.num_cols; c++)
504 eigen_cdf_diff = eigen_ch - eigen_cl;
509 float64_t h_bak = eigen_h(eigen_h.size()-1);
510 eigen_h(eigen_h.size()-1) = 0;
513 eigen_l2_plus_s2 = (eigen_s2.replicate(1,eigen_l.rows()).array().transpose().colwise() + eigen_l.array().pow(2)).matrix();
515 eigen_h2_plus_s2 = (eigen_s2.replicate(1,eigen_h.rows()).array().transpose().colwise() + eigen_h.array().pow(2)).matrix();
517 eigen_weighted_pdf_diff = (eigen_pl.array().colwise() * eigen_l.array() - eigen_ph.array().colwise() * eigen_h.array()).matrix();
520 eigen_h(eigen_h.size()-1) = h_bak;
virtual bool set_variational_distribution(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab)
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_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()
T * get_column_vector(index_t col) const
static float64_t univariate_log_pdf(float64_t sample, float64_t mu=0.0, float64_t sigma2=1.0)
virtual const char * get_name() const
all of classes and functions are contained in the shogun namespace
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 bool set_variational_distribution(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab)
virtual ~CLogitVGPiecewiseBoundLikelihood()
SGVector< float64_t > m_mu
virtual void init_likelihood()
SGVector< float64_t > m_s2