Skip to content

Commit 81c6511

Browse files
vandenmanvandenman
committed
can be tested
1 parent f719bdd commit 81c6511

File tree

2 files changed

+41
-22
lines changed

2 files changed

+41
-22
lines changed

‎src/ggm_model.cpp‎

Lines changed: 39 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -94,14 +94,6 @@ double GGMModel::log_density_impl_edge(size_t i, size_t j) const {
9494

9595
double trace_prod = -2 * (suf_stat_(j, j) * Uj2 + suf_stat_(i, j) * Ui2);
9696

97-
// This function uses the fact that the determinant doesn't change during edge updates.
98-
// double trace_prod = 0.0;
99-
// // TODO: we only need one of the two lines below, but it's not entirely clear which one
100-
// trace_prod += suf_stat_(j, j) * (omega_prop(j, j) - omega(j, j));
101-
// trace_prod += suf_stat_(i, i) * (omega_prop(i, i) - omega(i, i));
102-
// trace_prod += 2 * suf_stat_(i, j) * (omega_prop(i, j) - omega(i, j));
103-
// trace_prod - sum((aOmega_prop - aOmega) * SufStat)
104-
10597
double log_likelihood_ratio = (n_ * logdet - trace_prod) / 2;
10698
return log_likelihood_ratio;
10799

@@ -115,7 +107,7 @@ double GGMModel::log_density_impl_diag(size_t j) const {
115107
double cc12 = 1 - inv_omega_(j, j) * Uj2;
116108
double cc22 = 0 + Uj2 * Uj2 * inv_omega_(j, j);
117109

118-
double logdet = log(abs(cc11 * cc22 - cc12 * cc12));
110+
double logdet = std::log(std::abs(cc11 * cc22 - cc12 * cc12));
119111
double trace_prod = -2 * suf_stat_(j, j) * Uj2;
120112

121113
// This function uses the fact that the determinant doesn't change during edge updates.
@@ -182,9 +174,15 @@ void GGMModel::update_edge_parameter(size_t i, size_t j) {
182174
// double ln_alpha = log_density(omega_prop_) - log_density();
183175
double ln_alpha = log_density_impl_edge(i, j);
184176

185-
if (std::abs(ln_alpha - (log_density(omega_prop_) - log_density())) > 1e-6) {
186-
Rcpp::Rcout << "Warning: log density implementations do not match for edge (" << i << ", " << j << ")" << std::endl;
187-
}
177+
// {
178+
// double ln_alpha_ref = log_density(omega_prop_) - log_density();
179+
// if (std::abs(ln_alpha - ln_alpha_ref) > 1e-6) {
180+
// Rcpp::Rcout << "Warning: log density implementations do not match for edge (" << i << ", " << j << ")" << std::endl;
181+
// omega_.print(Rcpp::Rcout, "Current omega:");
182+
// omega_prop_.print(Rcpp::Rcout, "Proposed omega:");
183+
// Rcpp::Rcout << "ln_alpha: " << ln_alpha << ", ln_alpha_ref: " << ln_alpha_ref << std::endl;
184+
// }
185+
// }
188186

189187
ln_alpha += R::dcauchy(omega_prop_(i, j), 0.0, 2.5, true);
190188
ln_alpha -= R::dcauchy(omega_(i, j), 0.0, 2.5, true);
@@ -262,9 +260,16 @@ void GGMModel::update_diagonal_parameter(size_t i) {
262260
// 5) Acceptance ratio
263261
// double ln_alpha = log_density(omega_prop_) - log_density();
264262
double ln_alpha = log_density_impl_diag(i);
265-
if (std::abs(ln_alpha - (log_density(omega_prop_) - log_density())) > 1e-6) {
266-
Rcpp::Rcout << "Warning: log density implementations do not match for diag (" << i << ", " << i << ")" << std::endl;
267-
}
263+
// {
264+
// double ln_alpha_ref = log_density(omega_prop_) - log_density();
265+
// if (std::abs(ln_alpha - ln_alpha_ref) > 1e-6) {
266+
// Rcpp::Rcout << "Warning: log density implementations do not match for diag (" << i << ", " << i << ")" << std::endl;
267+
// // omega_.print(Rcpp::Rcout, "Current omega:");
268+
// // omega_prop_.print(Rcpp::Rcout, "Proposed omega:");
269+
// Rcpp::Rcout << "ln_alpha: " << ln_alpha << ", ln_alpha_ref: " << ln_alpha_ref << std::endl;
270+
// Rcpp::Rcout << "1e4 * diff: " << 10000 * (ln_alpha - ln_alpha_ref) << std::endl;
271+
// }
272+
// }
268273

269274
ln_alpha += R::dgamma(exp(theta_prop), 1.0, 1.0, true);
270275
ln_alpha -= R::dgamma(exp(theta_curr), 1.0, 1.0, true);
@@ -315,9 +320,16 @@ void GGMModel::update_edge_indicator_parameter_pair(size_t i, size_t j) {
315320

316321
// double ln_alpha = log_density(omega_prop_) - log_density();
317322
double ln_alpha = log_density_impl_edge(i, j);
318-
if (std::abs(ln_alpha - (log_density(omega_prop_) - log_density())) > 1e-6) {
319-
Rcpp::Rcout << "Warning: log density indicator implementations do not match for edge (" << i << ", " << j << ")" << std::endl;
320-
}
323+
// {
324+
// double ln_alpha_ref = log_density(omega_prop_) - log_density();
325+
// if (std::abs(ln_alpha - ln_alpha_ref) > 1e-6) {
326+
// Rcpp::Rcout << "Warning: log density implementations do not match for edge indicator (" << i << ", " << j << ")" << std::endl;
327+
// omega_.print(Rcpp::Rcout, "Current omega:");
328+
// omega_prop_.print(Rcpp::Rcout, "Proposed omega:");
329+
// Rcpp::Rcout << "ln_alpha: " << ln_alpha << ", ln_alpha_ref: " << ln_alpha_ref << std::endl;
330+
// }
331+
// }
332+
321333

322334
ln_alpha += std::log(1.0 - prior_inclusion_prob_(i, j)) - std::log(prior_inclusion_prob_(i, j));
323335

@@ -378,9 +390,15 @@ void GGMModel::update_edge_indicator_parameter_pair(size_t i, size_t j) {
378390

379391
// double ln_alpha = log_density(omega_prop_) - log_density();
380392
double ln_alpha = log_density_impl_edge(i, j);
381-
if (std::abs(ln_alpha - (log_density(omega_prop_) - log_density())) > 1e-6) {
382-
Rcpp::Rcout << "Warning: log density indicator implementations do not match for edge (" << i << ", " << j << ")" << std::endl;
383-
}
393+
// {
394+
// double ln_alpha_ref = log_density(omega_prop_) - log_density();
395+
// if (std::abs(ln_alpha - ln_alpha_ref) > 1e-6) {
396+
// Rcpp::Rcout << "Warning: log density implementations do not match for edge indicator (" << i << ", " << j << ")" << std::endl;
397+
// omega_.print(Rcpp::Rcout, "Current omega:");
398+
// omega_prop_.print(Rcpp::Rcout, "Proposed omega:");
399+
// Rcpp::Rcout << "ln_alpha: " << ln_alpha << ", ln_alpha_ref: " << ln_alpha_ref << std::endl;
400+
// }
401+
// }
384402
ln_alpha += std::log(prior_inclusion_prob_(i, j)) - std::log(1.0 - prior_inclusion_prob_(i, j));
385403

386404
// Prior change: add slab (Cauchy prior)

‎test_ggm.R‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
library(bgms)
22

33
# Dimension and true precision
4-
p <- 50
4+
p <- 10
55

66
adj <- matrix(0, nrow = p, ncol = p)
77
adj[lower.tri(adj)] <- rbinom(p * (p - 1) / 2, size = 1, prob = 0.3)
@@ -34,6 +34,7 @@ sampling_results <- bgms:::sample_ggm(
3434

3535
true_values <- zapsmall(Omega[upper.tri(Omega, TRUE)])
3636
posterior_means <- rowMeans(sampling_results[[2]]$samples)
37+
cbind(true_values, posterior_means)
3738

3839
plot(true_values, posterior_means)
3940
abline(0, 1)

0 commit comments

Comments
 (0)