Jan-28-2024, 08:53 AM
Hey, I have a problem with my homework. The goal of the first half is to sample 50 values of r_0,i from a Gaussian distribution with mean mu = 1 and sigma = 0.05. Also compute the 50 orbits corresponding to the initial conditions r(t_0) = r_0,i and dr/dt = 0, always considering L = 1.2. Finally compute the energy of each orbit, check that indeed it is conserved in time, and store its value E_i in an appropriate variable.
I think I've got the first half of the homework right:
I thought that everything is right but I always get the same error message:
I don't understand what I have to do, get my emcee.EnsembleSampler running :(
Thank you in advance.
I think I've got the first half of the homework right:
import numpy as np
from numpy import sqrt,exp,pi
import matplotlib.pyplot as plt
from scipy.io import readsav
from scipy.integrate import solve_ivp
import emcee
import corner
# 1.
A= 1
L0 = 1.2
r0 = 1.0
v0 = 0
psi0 = 0.0
t_span = (0, 50)
def gauss(L):
mu = 1.
sigma = 0.05
gaussian = exp(-0.5*((L-mu)/sigma)**2)/sqrt(2.*pi*sigma**2)
return gaussian
def kepler(t, a):
psi = L/a[0]**2
r = a[1]
v = (L**2/a[0]**3) - (A/a[0]**2)
return [r, v, psi]
a0 = (r0, v0, psi0)
def inverse_transform_sampling(n_samples=50):
L_values = np.linspace(0, 2, 1000)
CDF_values = np.cumsum(gauss(L_values)) / np.sum(gauss(L_values))
u_samples = np.random.uniform(0, 1, n_samples)
L_samples = np.interp(u_samples, CDF_values, L_values)
return L_samples
values = inverse_transform_sampling()
# 2. & 3
Ei = []
plt.figure(figsize=(10, 6))
for L in values:
sol = solve_ivp(kepler, t_span, a0, method='BDF')
x = sol.y[0]*np.cos(sol.y[2])
y = sol.y[0]*np.sin(sol.y[2])
energy = 0.5*sol.y[1]**2*0.5*(L/sol.y[0])**2-1/sol.y[0]
Ei.append(energy)
energy_change = np.max(np.abs(energy - Ei[-1]))
assert energy_change < 1e-10, f"Energy not conserved for orbit with L = {L}"
plt.plot(x, y, alpha=0.5)
plt.title('50 Orbits of Stars in Kepler Potential with Sampled Angular Momenta')
plt.xlim(-1.75, 1.5)
plt.ylim(-1.5, 1.5)
plt.show()
print(Ei)But now my problem comes. The second half of the homework we have to use emcee to fit a quadratic polinomial to the values of energy Ei as a function of the values of r_0,i. The error for each of our data points is E_err = 0.001.I thought that everything is right but I always get the same error message:
r0_samples = np.random.normal(1, 0.05, 50)
Ei_samples = np.array(Ei)
def model(theta, r):
a, b, c = theta
return a * r**2 + b * r + c
def log_prior(theta):
a, b, c = theta
if -1000 < a < 1000 and -1000 < b < 1000 and -1000 < c < 1000:
return 0.0
return -np.inf
def log_likelihood(theta, r, E, err):
model_vals = model(theta, r)
residual = E - model_vals
loglike = -0.5*np.sum(residual**2 / err**2 + np.log(err**2))
return loglike
def log_probability(theta, r, E, err):
lp = log_prior(theta)
if not np.isfinite(lp):
return -np.inf
logprob = lp + log_likelihood(theta, r, E, err)
return logprob
pos = [-2.4, 1.8, 3.6] + 1e-4*np.random.randn(32, 3)
nwalkers, ndim = pos.shape
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_probability, args=(r0_samples, Ei_samples, 0.001))
sampler.run_mcmc(pos, nsteps, progress=True)Error:ValueError Traceback (most recent call last)
Cell In[77], line 2
1 r0_samples = np.random.normal(1, 0.05, 50)
----> 2 Ei_samples = np.array(Ei)
4 def model(theta, r):
5 a, b, c = theta
ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (50,) + inhomogeneous part.I don't understand what I have to do, get my emcee.EnsembleSampler running :(
Thank you in advance.