thejoker icon indicating copy to clipboard operation
thejoker copied to clipboard

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

Open adrn opened this issue 3 years ago • 1 comments

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 avatar Jul 28 '20 13:07 adrn