I am trying to use emcee MCMC package to find alpha and characteristic mass of Schechter function. When I run the code, it always returns me the side boundary values from log_prior. Why is that? I add my code below and cannot find any solution... I have only observational data, which is defined as M (mass).
M = mass
# Define the Schechter mass function
def schechter_mass_function(M, alpha, Mc):
return M**(-alpha) * np.exp(-M / Mc)
# Define the log-likelihood function
def log_likelihood(theta, M):
alpha, Mc = theta
model = schechter_mass_function(M, alpha, Mc)
# Assuming Gaussian errors in your data
sigma = 0.1 # Adjust this based on your data
chi = -0.5 * np.sum((model - M)**2 / sigma**2)
return chi
# Define priors for alpha and Mc
def log_prior(theta):
alpha, Mc = theta
if (1.0 < alpha < 3.0) and (2.0 < Mc < 7.0):
return 0.0
return -np.inf
# Define the log-posterior function
def log_posterior(theta, M):
prior = log_prior(theta)
if not np.isfinite(prior):
return -np.inf
p = prior + log_likelihood(theta, M)
return p
# Set up emcee
ndim = 2 # Number of parameters
nwalkers = 100 # Number of walkers
nsteps = 200 # Number of steps for each walker
# Initialize walkers
initial_guess = np.array([2, 5]) + 0.1 * np.random.randn(nwalkers, ndim)
# Create the emcee sampler
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=(M,))
# Run the sampler
sampler.run_mcmc(initial_guess, nsteps, progress=True)
# Extract the samples
samples = sampler.chain[:, 100:, :].reshape((-1, ndim))
# Calculate the median and credible intervals for alpha and Mc
alpha_median, Mc_median = np.median(samples, axis=0)
alpha_credible_interval = np.percentile(samples[:, 0], [16, 84])
Mc_credible_interval = np.percentile(samples[:, 1], [16, 84])
I tried changing initial values, log_prior intervals, nwalkers, nsteps. Nothing helped to find the problem.