Closed ngierty closed 1 year ago
Could you post the full traceback?
Here's a screenshot of the traceback I get.
This one is really stumping me. Can you run the other tutorials ok?
Sorry for the delay. I have tried the Custom Model Tutorial but get the same issue on Spyder.
I'm able to somewhat get around the issue by following the example in the published paper and performing the fit in the command line. However, it is unclear to me how I would implement a custom likelihood if I can't perform the MCMC fit inside of Spyder.
I think this may have been an install issue on my end somewhere. In between doing a new installation of anaconda and reinstalling RadVel the issue cleared up. If I come across it again, I'll reopen the issue with hopefully some more insight.
I adapted the code in the SyntheticData.py to run in Spyder. However, the radvel.mcmc step fails with "IndexError: list index out of range". I'm running Python 3.9.7 with Spyder 5.4.2 on Mac OSX Ventura (13.0).
Here's the (slightly) modified code to reproduce the issue:
`''' Synthetic data tutorial for RadVel '''
import radvel import radvel.likelihood import numpy as np import matplotlib.pyplot as plt import copy from scipy import optimize
Initialize Keplerian Model Parameters
params = radvel.Parameters(1,basis='per tc secosw sesinw k') params['k1'] = radvel.Parameter(value=1) params['per1'] = radvel.Parameter(value=1) params['secosw1'] = radvel.Parameter(value=0.1) params['sesinw1'] = radvel.Parameter(value=+0.2) params['tc1'] = radvel.Parameter(value=0.) rv_mod = radvel.RVModel(params)
Generate Synthetic Data
t = np.random.random_sample(40) t = t * 4 - 2 t = np.sort(t) ti = np.linspace(-2,2,400) errvel = 0.3 jitter = 0.3 syn_rv = rv_mod(t) + np.random.normal(scale = np.sqrt(errvel2 + jitter2), size = t.size) plt.errorbar(t,syn_rv,yerr=errvel,fmt='.',label='Synthetic Data') plt.plot(ti,rv_mod(ti),label='Underlying Model') plt.xlabel('Time') plt.ylabel('RV') plt.legend() plt.show()
Generate the likelihood
like_syn = radvel.likelihood.RVLikelihood(rv_mod,t,syn_rv,np.zeros(t.size)+errvel) like_syn.params['gamma'] = radvel.Parameter(value=0) like_syn.params['jit'] = radvel.Parameter(value=errvel)
truths = copy.deepcopy(like_syn.params) # Store away model parameters for later reference like_syn.params.update(dict(k1=radvel.Parameter(value=3), secosw1=radvel.Parameter(value=0.1), sesinw1=radvel.Parameter(value=0.1), tc1=radvel.Parameter(value=0.1))) # perturb the starting guess
like_syn.params['jit'].vary = False # Don't vary jitter
Perform a maximum likelihood fit
res = optimize.minimize(like_syn.neglogprob_array, like_syn.get_vary_params(), method='Nelder-Mead' )
res = optimize.minimize(like_syn.neglogprob_array, like_syn.get_vary_params(), method='L-BFGS-B' )
print(res) print(like_syn) plt.errorbar(like_syn.x, like_syn.model(t)+like_syn.residuals(), yerr=like_syn.yerr, fmt='o') plt.plot(ti, like_syn.model(ti)) plt.xlabel('Time') plt.ylabel('RV') plt.show()
Instantiate a posterior object
post = radvel.posterior.Posterior(like_syn) post.params['per1'] = radvel.Parameter(value=1) post.params['k1'] = radvel.Parameter(value=1) post.params['jit'].vary = True post.priors += [radvel.prior.EccentricityPrior( 1 )] post.priors += [radvel.prior.Gaussian( 'jit', errvel, 0.1)] post.priors += [radvel.prior.Gaussian( 'per1', 1, 0.1)] post.priors += [radvel.prior.Gaussian( 'tc1', 0, 0.1)]
post.priors += [radvel.prior.SecondaryEclipsePrior(1, 0.5, 0.01)]
print(post)
Perform a maximum likelihood fit on posterior
res = optimize.minimize( post.neglogprob_array, post.get_vary_params(), method='Powell', options=dict(maxiter=100000,maxfev=100000,xtol=1e-8) )
plt.errorbar(post.likelihood.x, post.likelihood.model(t)+like_syn.residuals(), yerr=post.likelihood.yerr, fmt='o') plt.plot(ti, post.likelihood.model(ti)) plt.xlabel('Time') plt.ylabel('RV') plt.show()
Use MCMC to sample the posterior
df = radvel.mcmc(post, nrun=200, savename='rawchains.h5')`