谢克特功能司仪

Question

我正在尝试使用 emcee MCMC 包来查找 Schechter 函数的 alpha 和特征质量。 当我运行代码时，它总是返回 log_prior 的侧边界值。这是为什么？我在下面添加我的代码，但找不到任何解决方案... 我只有观测数据，其定义为 M（质量）。

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])

我尝试更改初始值、log_prior 间隔、nwalkers、nsteps。没有任何帮助找到问题。

Answer 1

我注意到在你对chi的定义中你有模型-M，但从数据来看它应该是模型-Phi（密度）。 它也可能有助于在 log(Phi) 空间中工作。