如何基于函数变量子集执行MCMC拟合？代码报错求助

我有一个含多参数的函数，示例为f(a,b,c)=a + b*x + c*x²，其中参数a、b、c可设为固定值或待拟合变量。我能正常同时拟合全部三个参数，但需要固定部分参数（如固定b拟合a和c，或固定b、c仅拟合a），原因是物理系统中部分参数已知，或数据量不足无法拟合过多参数。我已实现全参数拟合的可行代码，尝试用字典定义待拟合参数时，运行出现错误：ValueError: incompatible input dimensions (1, 50)，现寻求该问题的解决方法。

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)

问题出在**emcee不支持字典格式的参数输入**，它要求待拟合参数必须是**一维数组（ndarray）**格式，而你用字典初始化walker位置，导致维度不匹配报错。

要实现固定部分参数、拟合子集的需求，不需要用字典，而是通过以下思路修改代码：

1. 明确区分固定参数和待拟合参数，只把待拟合参数传入emcee的采样器
2. 在模型函数、似然、先验中，将固定参数和拟合参数合并使用

以下是修改后的可运行代码，以固定b，拟合a和c为例：

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

# -------------------------- 核心修改：区分固定和拟合参数 --------------------------
# 固定参数
fixed_params = {'b': true_b}
# 待拟合参数列表（按顺序）
fit_param_names = ['a', 'c']
ndim = len(fit_param_names)

# 定义带固定参数的模型函数
def f(x, fit_params, fixed_params):
    # 把拟合参数按对应名称赋值
    param_dict = dict(zip(fit_param_names, fit_params))
    # 合并固定参数
    param_dict.update(fixed_params)
    return param_dict['a'] + param_dict['b'] * x + param_dict['c'] * x**2

# 定义对数似然
def log_likelihood(fit_params, x, y, y_err, fixed_params):
    model = f(x, fit_params, fixed_params)
    return -0.5 * np.sum((y - model)**2 / y_err**2)

# 定义对数先验（只针对待拟合参数）
def log_prior(fit_params):
    # 对a和c设置均匀先验
    if all(-10.0 < p < 10.0 for p in fit_params):
        return 0.0
    return -np.inf

# 定义对数后验
def log_posterior(fit_params, x, y, y_err, fixed_params):
    lp = log_prior(fit_params)
    if not np.isfinite(lp):
        return -np.inf
    return lp + log_likelihood(fit_params, x, y, y_err, fixed_params)
# -----------------------------------------------------------------------------

# 初始化walker位置：只针对待拟合参数
nwalkers = 50
per = 0.1
# 待拟合参数的初始值
initial_fit_vals = [true_a, true_c]
initial_pos = [initial_fit_vals + per * np.random.randn(ndim) for _ in range(nwalkers)]

# 设置采样器：传入固定参数作为额外参数
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, 
                                args=(x_data, y_data, 0.05, fixed_params))

# 运行采样
nsteps = 2000
burnin = nsteps // 3
sampler.run_mcmc(initial_pos, nsteps, progress=True)

# 获取采样结果
samples = sampler.get_chain(discard=burnin, flat=True)

# 查看拟合结果：对应fit_param_names的顺序
print(f"拟合a的均值: {np.mean(samples[:,0]):.3f}")
print(f"拟合c的均值: {np.mean(samples[:,1]):.3f}")

灵活调整拟合参数的方法 #

如果需要更换拟合的参数子集，只需要修改两处：

• 修改fixed_params：把要固定的参数和对应值放入字典
• 修改fit_param_names：列出待拟合的参数名称，同时更新initial_fit_vals为对应初始值

比如要固定b和c，只拟合a，只需修改：

fixed_params = {'b': true_b, 'c': true_c}
fit_param_names = ['a']
initial_fit_vals = [true_a]

这种方式既符合emcee的参数要求，又能灵活实现固定部分参数的需求，避免了字典格式带来的维度问题。

内容的提问来源于stack exchange，提问作者Santiago del Palacio