adrn / thejoker

A custom Monte Carlo sampler for the (gravitational) two-body problem
MIT License
30 stars 8 forks source link

Can I use the same prior samples cache with different priors on linear parameters? #106

Open adrn opened 4 years ago

adrn commented 4 years ago

As originally asked by @astrosong:

If the sigma_K0 and sigma_v (and P_min, P_max) are different for each object, which means a different prior for each object, then can we still use the same prior_samples?

@adrn: Yes. The prior cache only contains the non-linear parameters (P, e, omega, M0, and optionally the extra variance parameter), so you can use different K and v priors with the same prior cache.

@astrosong: > @AstroSong

The second step: if the sigma_K0 and sigma_v (and P_min, P_max) are different for each object, which means a different prior for each object, then can we still use the same prior_samples?

Yes. The prior cache only contains the non-linear parameters (P, e, omega, M0, and optionally the extra variance parameter), so you can use different K and v priors with the same prior cache.

@adrn Sorry I am still a little confused. Do you mean that I use one prior to generate the shared prior cache, and input another prior to TheJoker for each star? The latter prior includes the parameter set same for the prior cache and the values for each star (sigma_K0, sigma_v0, P0, and v0_offsets).

Like Follows?

#Run once for all stars
with pm.Model():
    s = xu.with_unit(pm.Lognormal('s', 0, 0.5), u.m/u.s)
    e = xu.with_unit(pm.Deterministic('e', tt.constant(0)),
                              u.one)
    omega = xu.with_unit(pm.Deterministic('omega', tt.constant(0)),
                                       u.radian)
    prior = tj.JokerPrior.default(P_min=0.1*u.day, P_max=1000*u.day, 
                                                s=s, pars={'e': e, 'omega': omega})

prior_cache_file = 'prior_samples.hdf5'
if not os.path.exists(prior_cache_file):
    prior_samples = prior.sample(2*10**8)
    prior_samples.write(prior_cache_file)
#Run for each star
with pm.Model():
    dv0_1 = xu.with_unit(pm.Normal('dv0_1', 0, 5.), u.m/u.s)
    s = xu.with_unit(pm.Lognormal('s', 0, 0.5), u.m/u.s)
    e = xu.with_unit(pm.Deterministic('e', tt.constant(0)),
                     u.one)
    omega = xu.with_unit(pm.Deterministic('omega', tt.constant(0)),
                         u.radian)
    prior2 = tj.JokerPrior.default(P_min=0.1*u.day, P_max=1000*u.day, P0=1*u.day,
                                                v0_offsets=[dv0_1], s=s, 
                                                sigma_K0=30*u.km/u.s, sigma_v = 100*u.km/u.s,
                                                pars={'e': e, 'omega': omega})

prior_cache_file = 'prior_samples.hdf5'

with schwimmbad.MultiPool() as pool:
    joker = tj.TheJoker(prior2, pool=pool)
    samples = joker.rejection_sample(data, prior_cache_file )
adrn commented 4 years ago

@astrosong: Yes, what you have should almost work. But a few notes: