thejoker
thejoker copied to clipboard
Can I use the same prior samples cache with different priors on linear parameters?
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 )