adrn
commented
2 years ago Hi @AstroJLin - here's a small example: https://gist.github.com/56707a370423d883a33bcc9049bdbf7a
Open AstroJLin opened 2 years ago
adrn
commented
2 years ago Hi @AstroJLin - here's a small example: https://gist.github.com/56707a370423d883a33bcc9049bdbf7a
AstroJLin
commented
2 years ago Hi, Here, I fix the period with P = xu.with_unit(pm.Constant('P', P), u.day). However, it does not seen to run and the error is as follows. Should I how to deal with this error?
 File "/home/jlin/xinglong/RV/fixP.py", line 110, in
builtins.ValueError: cannot reshape array of size 0 into shape (0,newaxis)
Thanks!
adrn
commented
2 years ago What versions of cython, numpy, and thejoker do you have installed?
AstroJLin
commented
2 years ago Hi, Here, I use cython (0.29.21), numpy(1.21.1) and thejoker (1.1). In fact, If I run a small example that you give me, it can work. However, when I use my data to run with P = xu.with_unit(pm.Constant('P', P), u.day), this error will appear. So, I do not know what is happen？
Thanks!
adrn
commented
2 years ago What type is your variable P? You should probably name it something else, like period or _P since you redefine it as a model parameter. What if you try:
P = xu.with_unit(pm.Constant('P', float(_P)), u.day)?
AstroJLin
commented
2 years ago Hi,
Although have named it something else, it still reports this error. Here I attach my data and code, could you help me to see what is happen?
Thanks!
from astropy.io import ascii from astropy.time import Time import astropy.units as u import matplotlib.pyplot as plt
import numpy as np
import pymc3 as pm import exoplanet.units as xu import exoplanet as xo from exoplanet.distributions import Angle import corner import arviz as az
import thejoker as tj import emcee seed = 42 rnd = np.random.default_rng(seed=seed)
tbl = ascii.read(
"""
2458183.36136016 -97.2 3.1
2458183.37733315 -37.7 4.9
2458183.39330614 24.5 3.8
2458187.34557366 -7.6 3
2458187.36154661 -89.7 2.5
2458187.37821403 -155.2 2
2458187.39418697 -193.7 2
2459295.33498398 141.1 2.2
2459295.35165123 153.2 2.1
2459295.36762401 147 1.9
2457091.33443744 -187.9 1.3
2458902.39203349 143 0.6
""",
names=['BJD', 'rv', 'rv_err'])
tbl['rv'].unit = u.km/u.s
tbl['rv_err'].unit = u.km/u.s
data = tj.RVData( t=Time(tbl['BJD'], format='jd', scale='tcb'), rv=u.Quantity(tbl['rv']), rv_err=u.Quantity(tbl['rv_err']))
data.plot(markersize=3) plt.show() period = 0.2556686 with pm.Model() as model: e = xu.with_unit(pm.Constant('e', 0),u.one) P = xu.with_unit(pm.Constant('P', float(period)), u.day) prior = tj.JokerPrior.default(
#P_min=0.1*u.day, P_max=1.*u.day,
sigma_K0=170*u.km/u.s, # scale of the prior on semiamplitude, K
sigma_v=50*u.km/u.s, # std dev of the prior on the systemic velocity, v0
#v0_offsets=[dv0_1],
pars={'e':e,'P':P}
)joker = tj.TheJoker(prior) samples = joker.rejection_sample(data,prior_samples=100_000, max_posteriorsamples=128) = tj.plot_rv_curves(samples, data=data)
adrn
commented
2 years ago Oh interesting - there is something weird about the way pm.Constant() is storing the value (it keeps it internally as an int, so the period becomes 0, which is invalid).
Try changing the definitions of e and P to:
e = xu.with_unit(pm.Deterministic('e', tt.as_tensor_variable(0)), u.one)
P = xu.with_unit(pm.Deterministic('P', tt.as_tensor_variable(period)), u.day)Also added your example to the bottom here: https://gist.github.com/56707a370423d883a33bcc9049bdbf7a
AstroJLin
commented
2 years ago Hi Adrian, Thank you very much for your help and reply.
I have been able to fixed the period after using your suggestions. And then, when I try to run the pmx.sample() in the mcmc, it throw an error (builtins.TypeError: ('Wrong number of dimensions: expected 0, got 1 with shape (256,).', 'Container name "K"')). I do not know what is causing this error. Could you give me some suggetions? Here I attach code and this error.
###
with pm.Model():
e = xu.with_unit(pm.Deterministic('e', tt.constant(0)),
P = xu.with_unit(pm.Deterministic('P', tt.constant(period)), u.day)
prior_mcmc = tj.JokerPrior.default(
sigma_K0=100*u.km/u.s,
sigma_v=50*u.km/u.s,
pars={'e': e, 'P':P}
)
joker_mcmc = tj.TheJoker(prior_mcmc, random_state=rnd)
mcmc_init = joker_mcmc.setup_mcmc(data, samples)
print (mcmc_init)
(an error??? ) trace = pmx.sample(
tune=500, draws=1000,
start=mcmc_init,
random_seed=seed,
cores=1, chains=2)
print (az.summary(trace, var_names=prior_mcmc.par_names))
mcmc_samples = joker_mcmc.trace_to_samples(trace, data=data)
mcmc_samples.wrap_K()
print (mcmc_samples.mean().tbl['P', 'e', 'omega', 'M0', 's', 'K', 'v0'])
df = mcmc_samples.tbl.to_pandas()
df = df.drop(columns=['e', 's', 'omega','M0'])
_ = corner.corner(df,quantiles=(0.16,0.5,0.84),show_titles=True, title_kwargs={"fontsize": 12})
plt.show()
fig = plt.figure(figsize=(10, 8))
_ = tj.plot_phase_fold(mcmc_samples.median(), data)
plt.xlabel('Phase')
plt.ylabel('Radial Velocity (km/s)')
plt.show()
###File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/theano/link/basic.py", line 82, in set self.storage[0] = self.type.filter_inplace( File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/theano/graph/type.py", line 95, in filter_inplace raise NotImplementedError()
builtins.NotImplementedError
During handling of the above exception, another exception occurred:
 File "/home/jlin/xinglong/RV/fixP.py", line 161, in
builtins.TypeError: ('Wrong number of dimensions: expected 0, got 1 with shape (256,).', 'Container name "K"')
adrn
commented
2 years ago For your data, thejoker returns many samples -- i.e., your samples variable has 256 posterior samples. When you run setup_mcmc(), you should only pass one sample. Here you can do that by just picking one at random, e.g.,
mcmc_init = joker_mcmc.setup_mcmc(data, samples[0])or by taking the MAP or maximum likelihood sample, e.g.,
loglike = samples.ln_unmarginalized_likelihood(data)
mcmc_init = joker_mcmc.setup_mcmc(data, samples[loglike.argmax()])Hi Adrian, Thank you very much for your reply again.
Here, I have two questions.
First, when I using the maximum likelihood sample, I do not use the median MCMC sample to fold the data and plot residuals relative to our inferred RV model. Here, an error is report as follow:
 File "/home/jlin/xinglong/RV/fixP.py", line 202, in
which means that I should set up the number of samples ( mcmc_samples.median()[?]). However, I do not know which samples corresponding to the median MCMC sample?
Second, If I do not try to fix the period, I do not need to set up the maximum likelihood sample as an example from the Thompson et al. 2019. So, I want to know the reason for this.
Thanks!
Hi,
I want to fix the period, which should how to set up?
Thanks!