1

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

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

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

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

4

3 回答 3

3

由于您使用的是标准优化例程，因此您可以将要最小化的函数表示为对数似然。如果它表示为对数似然，并且您要做的只是探索参数的后验分布，您可能会发现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() 插入司​​仪：

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)
于 2014-12-09T23:13:58.527 回答
2

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

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

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

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

（这都是 PyMC 2.3）

于 2014-11-30T23:47:20.770 回答
0

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

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

于 2014-12-31T00:00:20.867 回答