emcee for Schechter function

32 Views Asked by At

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.

0

There are 0 best solutions below