如何仅使用函数变量的子集进行 MCMC 拟合？

Question

我有一个有很多参数/变量的函数。为简单起见，我们可以假设 f(a,b,c)=a + b * x + c * x**2，其中参数“a”、“b”和“c”可以是固定参数或变量我想要适合。 我可以毫无问题地同时将 f(a,b,c) 与三个参数拟合。问题是当我想修复“b”并适合“a”和“c”，或者修复“b”和“c”并只适合“a”时。动机是，在物理系统中，其中一些参数是已知的，因此我可以专注于拟合其他参数（或者也许我没有足够的数据点来拟合许多参数，我更喜欢修复一些参数并只释放最相关的）。

这是一个有效的代码版本，修复以适应所有参数：

import numpy as np
import emcee
import corner
import matplotlib.pyplot as plt

# True parameters of the function
true_a = 1.0
true_b = 2.0
true_c = 0.5

# Generate synthetic data with noise
np.random.seed(42)
x_data = np.linspace(-8, 10, 50)
y_true = true_a + true_b * x_data + true_c * x_data**2
y_data = y_true + np.random.normal(scale=0.05, size=len(x_data))

# Define the quadratic function
def f(x, a, b, c):
    return a + b * x + c * x**2

# Define the log-likelihood function
def log_likelihood(params, x, y, y_err):
    a, b, c = params
    model = f(x, a, b, c)
    return -0.5 * np.sum((y - model)**2 / y_err**2)

# Define the log-prior function
def log_prior(params):
    # Uniform priors for a, b, c
    if all(-10.0 < p < 10.0 for p in params):
        return 0.0
    return -np.inf

# Define the log-posterior function
def log_posterior(params, x, y, y_err):
    lp = log_prior(params)
    if not np.isfinite(lp):
        return -np.inf
    return lp + log_likelihood(params, x, y, y_err)

# Perform MCMC fitting
ndim = 3  # Number of parameters (a, b, c)
nwalkers = 50
nsteps = 2000
burnin = nsteps // 3  # Burn-in steps (1/3 of the total steps)

# Initialize walkers with random values around the true parameters
initial_params = [true_a, true_b, true_c]
per = 0.1
initial_pos = [initial_params + per * np.random.randn(ndim) for _ in range(nwalkers)]

# Set up the MCMC sampler
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=(x_data, y_data, 0.5))

# Run the sampler with burn-in
sampler.run_mcmc(initial_pos, nsteps, progress=True)

# Get the samples after burn-in
samples = sampler.get_chain(discard=burnin, flat=True)

现在这是我尝试定义哪些参数适合使用字典作为参数（尽管这可能不是正确的方法？）

import numpy as np
import emcee
import corner
import matplotlib.pyplot as plt

# True parameters of the function
true_a = 1.0
true_b = 2.0
true_c = 0.5

# Generate synthetic data with noise
np.random.seed(42)
x_data = np.linspace(-8, 10, 50)
y_true = true_a + true_b * x_data + true_c * x_data**2
y_data = y_true + np.random.normal(scale=0.05, size=len(x_data))

# Define the function
def f(x, params):
    a = params.get('a', 0.0)
    b = params.get('b', 0.0)
    c = params.get('c', 0.0)
    return a + b * x + c * x**2

# Define the log-likelihood function
def log_likelihood(params, x, y, y_err):
    model = f(x, params)
    return -0.5 * np.sum((y - model)**2 / y_err**2)

# Define the log-prior function
def log_prior(params):
    # Uniform priors for a, b, c
    a = params.get('a', 0.0)
    b = params.get('b', 0.0)
    c = params.get('c', 0.0)
    if all(-10.0 < p < 10.0 for p in [a, b, c]):
        return 0.0
    return -np.inf

# Define the log-posterior function
def log_posterior(params, x, y, y_err):
    lp = log_prior(params)
    if not np.isfinite(lp):
        return -np.inf
    return lp + log_likelihood(params, x, y, y_err)


# Initialize walkers with random values around the true parameters
nwalkers = 50
initial_params = {'a': true_a, 'b': true_b, 'c': true_c} # here I would like to have any combination of parameters, not necessarily all three
per = 0.1
initial_pos = [{key: val + per * np.random.randn() for key, val in initial_params.items()} for _ in range(nwalkers)]

# Perform MCMC fitting
ndim = len(initial_params)  # Number of parameters 
nsteps = 2000
burnin = nsteps // 3  # Burn-in steps (1/3 of the total steps)


# Set up the MCMC sampler
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=(x_data, y_data, 0.05))

# Run the sampler with burn-in
sampler.run_mcmc(initial_pos, nsteps, progress=True)

# Get the samples after burn-in
samples = sampler.get_chain(discard=burnin, flat=True)

但是当它运行时它会说： ---> 63 Sampler.run_mcmc(initial_pos, nsteps, Progress=True) ValueError：输入尺寸不兼容 (1, 50)

我非常感谢有关如何解决此问题的反馈。

Answer 1

实现此目的的一种方法是使用 bilby 包，它为 emcee 提供了一个包装器。使用 bilby，您可以设置模型的任何固定参数以获得

DeltaFunction

先验。例如：

import numpy as np

from bilby import run_sampler
from bilby.core.likelihood import GaussianLikelihood
from bilby.core.prior import Uniform, DeltaFunction


# model function
def model(x, a, b, c):
    return a + b * x + c * x**2


# True parameters of the function
true_a = 1.0
true_b = 2.0
true_c = 0.5

# Generate synthetic data with noise
np.random.seed(42)
x_data = np.linspace(-8, 10, 50)
y_true = model(x_data, true_a, true_b, true_c)

y_err = 0.05
y_data = y_true + np.random.normal(scale=y_err, size=len(x_data))

# set the Gaussian likelihood
likelihood = GaussianLikelihood(x_data, y_data, model, sigma=y_err)

# set the priors (uniform on a and b, but fixed for c)
prior = {
    "a": Uniform(-10, 10, "a", latex_label="$a$"),
    "b": Uniform(-10, 10, "b", latex_label="$b$"),
    "c": DeltaFunction(true_c, "c"),  # c is fixed at the true value
}

nwalkers = 50
nsteps = 2000
burnin = nsteps // 3

results = run_sampler(
    likelihood=likelihood,
    priors=prior,
    sampler="emcee",
    nwalkers=50,
    nsteps=2000,
    nburn=burnin,
    save=False,
)

# a pandas DataFrame containing the posterior samples
samples = results.posterior

# plot the posterior samples (for the non-fixed parameters) on a corner plot
results.plot_corner()

输出图（默认保存到名为

outdir

的目录中）为：