带有 pymc 的参数向量的非线性回归答案

【问题标题】：non-linear regression with a parameter vector with pymc带有 pymc 的参数向量的非线性回归
【发布时间】：2015-01-17 16:57:33
【问题描述】：

我目前使用 scipy.optimize.minimize 和 scipy.optimize.leastsq 对我的数据集执行非线性回归。我想使用 PyMC(3) 来研究拟合过程中涉及的所有参数的后验。我在 SO 上遇到了这个previous answer。

这是一个很好的例子，我看到的大多数其他例子都是线性回归的。但是，该示例并不完全适合我的目的。我的模型具有可变数量的参数，我将拟合其中的一个子集。该子集通常在 1 到 20 个参数的范围内，但有时更多。使用 scipy 最小化器，这些不同的参数以 1D np.ndarray 的形式传递给成本函数，p，例如

def chi2(p, *args):
    xdata = args[0]
    return p[0] + xdata * p[1] + ........

在上面给出的链接中，@pymc.deterministic 装饰的高斯函数具有关键字参数。这对我来说是不切实际的，因为同一个代码块需要处理不同（和相当大）数量的参数。有没有办法提供参数向量？我还必须为每个参数提供一个先验列表。但是，我有每个参数 [(min, max)...] 的下限和上限列表，所以这不是问题。

【问题讨论】：

参数的数量不应该随着迭代而改变，不是吗？你能详细说明一下吗？
克里斯，我正在尝试将 pymc 固定到 lmfit-py 项目中。这个想法是使用 pymc 对曲线拟合场景中的参数执行贝叶斯分析。在 lmfit 中，创建了一个返回曲线拟合残差的目标函数。从那以后，我已经弄清楚了如何在这种情况下使用 pymc，这很棒。本质上，我创建了一个@pm.observed likelihood 函数，它根据用户目标函数返回一个pm.normal_like 分布。在 lmfit 中，您的曲线拟合可能有 M 个参数，您允许其中 N 个参数变化，0

标签： pymc pymc3

【解决方案1】：

对于一组参数，通常最好使用向量并取点积：

@deterministic
def theta(beta=beta):
    return beta.dot(X)

或者如果你想要一个明确的基线均值：

@deterministic
def theta(mu=mu, beta=beta):
    return mu + beta.dot(X)

（这都是 PyMC 2.3）

【讨论】：

【解决方案2】：

由于您使用的是标准优化例程，您可以将要最小化的函数表示为对数似然。如果它表示为对数似然，并且您只想探索参数的后验分布，您可能会发现 emcee 更易于使用。就我而言，在开始研究 mcmc 方法之前，我实际上是在最小化对数似然。

from scipy.optimize import minimize

bnds = ((.0001,20),(.0001,20),(.0001,20),(.0001,20))

solution = minimize(my_function_ll, np.array([1.1,1.2,1.3,1.4]), 
                    args=(my_data),
                    bounds=bnds )

我所要做的就是像这样将 my_function_ll() 插入 emcee：

import emcee

# for reproducible results
np.random.seed(0)

ndim = 4  # number of parameters in the model
nwalkers = 10  # number of MCMC walkers
nburn = 1000  # "burn-in" period to let chains stabilize
nsteps = 10000  # number of MCMC steps to take

# we'll start at random locations between 0 and 2
starting_guesses = 2*np.random.rand(nwalkers, ndim)

def log_prior(theta):
    x1,x2,x3,x4 = theta
    # here you can specify boundary. In my case I was using 0 - 20
    if 0.0001 < x1 < 20 and 0.0001 < x2 < 20 and 0.0001 < x3 < 20 and 0.0001 < x4 < 20:
        return 0.0
    return -np.inf

def log_posterior(theta, observation_array):
    lp = log_prior(theta)
    if not np.isfinite(lp):
        return -np.inf    
    return log_prior(theta) + my_function_ll(theta, observation_array) 

sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, 
                                args=[my_data])
sampler.run_mcmc(starting_guesses, nsteps)

sample = sampler.chain[:, nburn:, :].reshape(-1, 4)

现在您可以比较 MCMC 和 MLE 结果：

print "MLE: ", solution.x
print "MCMC: ", sample.mean(0)

【讨论】：

【解决方案3】：

以下链接包含很多我可以用来开始的信息。

https://groups.google.com/forum/#!msg/pymc/WUg0tZjNAL8/jsVRFe1DMT4J

【讨论】：