Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
adrn opened this issue Jul 28, 2020 · 1 comment
Labels

Comments

@adrn
Copy link
Owner

adrn commented Jul 28, 2020

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
Copy link
Owner Author

adrn commented Jul 28, 2020

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

  • Admittedly, this is a bit klugy (having to redefine and pass in the same nonlinear parameters in each model). I should think about allowing specifying only the linear parameter priors to initialize JokerPrior. In the second pm.Model() block it doesn't even use the parameters you define because you rejection sample from a cache file, so I guess it should be possible to initialize a JokerPrior under the assumption that you will sample from a pre-cached file.
  • In your first pm.Model block, you still have to specify sigma_K0 and sigma_v. You can either set them to arbitrary values, or (since they aren't used to generate the prior cache) you could set them to nan like: sigma_v=np.nan*u.km/u.s, sigma_K0=np.nan*u.km/u.s.
  • The prior on dv0_1 is only needed if you are working with data that comes from two different sources, with a possible constant offset between the two sources (e.g., telescopes). If you have data from one source, you can ignore that parameter.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

1 participant