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

Issues with Parameter Extraction from Emcee #522

Open
AjlaTr opened this issue Jun 4, 2024 · 0 comments
Open

Issues with Parameter Extraction from Emcee #522

AjlaTr opened this issue Jun 4, 2024 · 0 comments

Comments

@AjlaTr
Copy link

AjlaTr commented Jun 4, 2024

My team and I were hoping to extract the inclination angle for an SB1 binary system using emcee. We ran it with 5 parameters: inclination, period, the semi-amplitude of radial velocity, period, eccentricity, periastron, and center-of-mass velocity. However, when we did so, we noticed that the value of inclination we received was too small compared to the actual value. Additionally, the corner plot has a large distribution of values.
When I attempted to increase the number of walkers, the value for inclination slightly improved, but this caused issues with the other 4 parameters. As I am new to emcee, I'm not certain how to proceed.

`def binary_model(params, t):
K, P, e, omega, i = params
M = 2 * np.pi * t / P
f = M + e * np.sin(M) # Kepler's equation approximation for eccentric anomaly
rv_pred = K * (np.cos(f + omega) + e * np.cos(omega)) * np.sin(np.radians(i))
return rv_pred

def log_likelihood(params, t, rv, rv_err):
rv_pred = binary_model(params, t)
sigma2 = rv_err**2
return -0.5 * np.sum((rv - rv_pred)**2 / sigma2 + np.log(sigma2))

def log_prior(params):
K, P, e, omega, i = params
if 0 < K < 100 and 0 < P < 100 and 0 <= e < 1 and 0 <= omega < 2*np.pi and 0 < i < 90:
return np.log(np.sin(np.radians(i))) # Prior proportional to sin(i)
return -np.inf

def log_posterior(params, t, rv, rv_err):
lp = log_prior(params)
if not np.isfinite(lp):
return -np.inf
return lp + log_likelihood(params, t, rv, rv_err)

nwalkers = 100
ndim = 5 # Number of parameters in the model
initial_params = [65.49682, 59.93311, 0.01499, 0.04211, 60] # Initial guesses for the parameters
p0 = [initial_params + 1e-4*np.random.randn(ndim) for i in range(nwalkers)]
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=(t, rv_obs, rv_err))

sampler.run_mcmc(p0, 10000) # Increase the number of steps

samples = sampler.get_chain(discard=2000, thin=15, flat=True) # Discard a larger initial portion
incl_samples = samples[:, 3] # Extract inclination samples`

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

No branches or pull requests

1 participant