|
1 | 1 | # pylint: skip-file |
2 | | -import chex |
| 2 | + |
| 3 | +import distrax |
| 4 | +import haiku as hk |
| 5 | +import optax |
| 6 | +from jax import numpy as jnp |
| 7 | +from surjectors import Chain, MaskedCoupling, TransformedDistribution |
| 8 | +from surjectors.conditioners import mlp_conditioner |
| 9 | +from surjectors.util import make_alternating_binary_mask |
| 10 | + |
| 11 | +from sbijax import SNL |
| 12 | + |
| 13 | + |
| 14 | +def prior_model_fns(): |
| 15 | + p = distrax.Independent( |
| 16 | + distrax.Uniform(jnp.full(2, -3.0), jnp.full(2, 3.0)), 1 |
| 17 | + ) |
| 18 | + return p.sample, p.log_prob |
| 19 | + |
| 20 | + |
| 21 | +def simulator_fn(seed, theta): |
| 22 | + p = distrax.MultivariateNormalDiag(theta, 0.1 * jnp.ones_like(theta)) |
| 23 | + y = p.sample(seed=seed) |
| 24 | + return y |
| 25 | + |
| 26 | + |
| 27 | +def log_density_fn(theta, y): |
| 28 | + prior = distrax.Uniform(jnp.full(2, -3.0), jnp.full(2, 3.0)) |
| 29 | + likelihood = distrax.MultivariateNormalDiag( |
| 30 | + theta, 0.1 * jnp.ones_like(theta) |
| 31 | + ) |
| 32 | + |
| 33 | + lp = jnp.sum(prior.log_prob(theta)) + jnp.sum(likelihood.log_prob(y)) |
| 34 | + return lp |
| 35 | + |
| 36 | + |
| 37 | +def make_model(dim): |
| 38 | + def _bijector_fn(params): |
| 39 | + means, log_scales = jnp.split(params, 2, -1) |
| 40 | + return distrax.ScalarAffine(means, jnp.exp(log_scales)) |
| 41 | + |
| 42 | + def _flow(method, **kwargs): |
| 43 | + layers = [] |
| 44 | + for i in range(2): |
| 45 | + mask = make_alternating_binary_mask(dim, i % 2 == 0) |
| 46 | + layer = MaskedCoupling( |
| 47 | + mask=mask, |
| 48 | + bijector=_bijector_fn, |
| 49 | + conditioner=mlp_conditioner([8, 8, dim * 2]), |
| 50 | + ) |
| 51 | + layers.append(layer) |
| 52 | + chain = Chain(layers) |
| 53 | + base_distribution = distrax.Independent( |
| 54 | + distrax.Normal(jnp.zeros(dim), jnp.ones(dim)), |
| 55 | + 1, |
| 56 | + ) |
| 57 | + td = TransformedDistribution(base_distribution, chain) |
| 58 | + return td(method, **kwargs) |
| 59 | + |
| 60 | + td = hk.transform(_flow) |
| 61 | + td = hk.without_apply_rng(td) |
| 62 | + return td |
3 | 63 |
|
4 | 64 |
|
5 | 65 | def test_snl(): |
6 | | - chex.assert_equal(1, 1) |
| 66 | + rng_seq = hk.PRNGSequence(0) |
| 67 | + y_observed = jnp.array([-1.0, 1.0]) |
| 68 | + |
| 69 | + prior_simulator_fn, prior_logdensity_fn = prior_model_fns() |
| 70 | + fns = (prior_simulator_fn, prior_logdensity_fn), simulator_fn |
| 71 | + |
| 72 | + snl = SNL(fns, make_model(2)) |
| 73 | + params, info = snl.fit( |
| 74 | + next(rng_seq), |
| 75 | + y_observed, |
| 76 | + n_rounds=1, |
| 77 | + optimizer=optax.adam(1e-4), |
| 78 | + sampler="slice", |
| 79 | + ) |
| 80 | + _ = snl.sample_posterior(params, 2, 100, 50, sampler="slice") |
0 commit comments