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.
Why do I get this error message?I don't understand what I have to do, get my emcee.EnsembleSampler running :(
Thank you in advance.