use fugue::inference::abc::{abc_smc_weighted, ABCSMCConfig, EuclideanDistance};
use fugue::inference::vi::{optimize_meanfield_vi_with_config, Support, VIConfig};
use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;
const POSTERIOR_MEAN: f64 = 1.2;
const POSTERIOR_VAR: f64 = 0.2;
const POSTERIOR_SD: f64 = 0.4472135954999579;
fn observed_model() -> Model<f64> {
sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap())
.bind(|mu| observe(addr!("y"), Normal::new(mu, 0.5).unwrap(), 1.5).map(move |_| mu))
}
fn forward_sim_model() -> Model<f64> {
sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap())
.bind(|mu| sample(addr!("y_sim"), Normal::new(mu, 0.5).unwrap()).map(move |_| mu))
}
#[test]
fn fg25_smc_example_recovers_known_posterior() {
let mut rng = StdRng::seed_from_u64(42);
let config = SMCConfig {
resampling_method: ResamplingMethod::Systematic,
ess_threshold: 0.5,
rejuvenation_steps: 3,
};
let result = adaptive_smc(&mut rng, 2000, observed_model, config);
let mean: f64 = result
.iter()
.filter_map(|p| p.trace.get_f64(&addr!("mu")).map(|mu| p.weight * mu))
.sum();
let var: f64 = result
.iter()
.filter_map(|p| {
p.trace
.get_f64(&addr!("mu"))
.map(|mu| p.weight * (mu - mean).powi(2))
})
.sum();
assert!(
(mean - POSTERIOR_MEAN).abs() < 0.15,
"SMC mean {mean} deviates from exact posterior mean {POSTERIOR_MEAN}"
);
assert!(
(var - POSTERIOR_VAR).abs() < 0.1,
"SMC var {var} deviates from exact posterior var {POSTERIOR_VAR}"
);
assert!(result.log_evidence.is_finite());
}
#[test]
fn fg25_abc_example_recovers_known_posterior() {
let observed: Vec<f64> = vec![1.5];
let mut rng = StdRng::seed_from_u64(7);
let config = ABCSMCConfig {
initial_tolerance: 2.0,
tolerance_schedule: vec![1.0, 0.5, 0.25, 0.1],
particles_per_round: 500,
};
let result = abc_smc_weighted(
&mut rng,
forward_sim_model,
|trace| vec![trace.get_f64(&addr!("y_sim")).unwrap()],
&observed,
&EuclideanDistance,
config,
200_000,
)
.expect("ABC-SMC should complete with this many particles/attempts");
let mean = result
.weighted_mean(&addr!("mu"))
.expect("mu present in every particle");
assert!(
(mean - POSTERIOR_MEAN).abs() < 0.3,
"ABC-SMC mean {mean} deviates from target posterior mean {POSTERIOR_MEAN}"
);
}
#[test]
fn fg25_vi_example_recovers_known_posterior() {
let mut rng = StdRng::seed_from_u64(11);
let mut guide = MeanFieldGuide::new();
guide.add_latent(addr!("mu"), Support::Real, 0.0);
let config = VIConfig {
n_iterations: 800,
n_samples_per_iter: 32,
base_learning_rate: 0.3,
..VIConfig::default()
};
let result = optimize_meanfield_vi_with_config(&mut rng, observed_model, guide, &config);
let VariationalParam::Normal { mu, log_sigma } = result
.guide
.params
.get(&addr!("mu"))
.expect("guide has a factor for mu")
else {
panic!("Support::Real latent must produce a Normal factor");
};
let sigma = log_sigma.exp();
assert!(
(*mu - POSTERIOR_MEAN).abs() < 0.15,
"VI mean {mu} deviates from exact posterior mean {POSTERIOR_MEAN}"
);
assert!(
(sigma - POSTERIOR_SD).abs() < 0.15,
"VI sd {sigma} deviates from exact posterior sd {POSTERIOR_SD}"
);
}