标签: pymc3

我试图用并行化加速pymc3采样,我看到只有适度的好处.

在i7 MacBook Pro上,我能够将总运行时间从25分钟(njobs = 1)减少到13分钟(njobs = 6).由于在pymc实际开始采样之前需要大约4分钟,因此增加相对较小.

问题是 - 是否有人成功地使用带有pymc3的GPU以及需要6-8分钟采样的模型可以获得多少好处？(我的MacBook有nvidia GT 750M 2Gb)

出于比较目的，我想利用 PyMC3 之外的后验密度函数。

对于我的研究项目，我想了解 PyMC3 与我自己定制的代码相比的性能如何。因此，我需要将它与我们自己的内部采样器和似然函数进行比较。

我想我想出了如何调用内部PyMC3后验的方法，但是感觉很别扭，想知道是否有更好的方法。现在我正在手动转换变量，而我应该能够将参数字典传递给 pymc 并获得后验密度。这可能以直接的方式进行吗？

非常感谢！

演示代码：

import numpy as np
import pymc3 as pm
import scipy.stats as st

# Simple data, with sigma = 4. We want to estimate sigma
sigma_inject = 4.0
data = np.random.randn(10) * sigma_inject

# Prior interval for sigma
a, b = 0.0, 20.0

# Build PyMC model
with pm.Model() as model:
    sigma = pm.Uniform('sigma', a, b)      # Prior uniform between 0.0 and 20.0
    likelihood = pm.Normal('data', 0.0, sd=sigma, observed=data)

# …

是否可以通过简单的方式从用户提供的 PyMC3 目标度量中进行采样？即我希望能够提供黑盒函数logposterior(theta)以及grad_logposterior(theta)那些和样本，而不是在 PyMC3s 建模语言中指定模型。

我想使用pymc3来估计霍奇金·赫x黎神经元模型中的未知参数和状态。我在pymc中的代码基于http://healthyalgorithms.com/2010/10/19/mcmc-in-python-how-to-stick-a-statistical-model-on-a-system-dynamics-model- in-pymc /并执行得很好。

#parameter priors
@deterministic
def HH(priors in here)
    #model equations
    #return numpy arrays that somehow contain the probability distributions as elements
    return V,n,m,h

#Make V deterministic in one line. Seems to be the magic that makes this work.
V = Lambda('V', lambda HH=HH: HH[0])

#set up the likelihood
A = Normal('A',mu=V,tau=2.0,value=voltage_data,observed = True)
#run the sampling...

在pymc3中，我无法使用Lambda技巧。这是我的尝试之一：

import numpy as np
import theano.tensor as tt
from pymc3 import Model, Normal, Uniform, Deterministic, sample, NUTS, Metropolis, find_MAP …

我一直在关注 PyMC3 的高斯混合模型示例：https : //github.com/pymc-devs/pymc3/blob/master/pymc3/examples/gaussian_mixture_model.ipynb 并让它与人工数据集很好地工作。

我已经用一个真实的数据集尝试过它，我正在努力让它给出合理的结果：

关于我应该缩小/扩大/改变哪些参数以获得更好的拟合的任何想法？痕迹似乎很稳定。这是我从示例中调整的模型片段：

model = pm.Model()
with model:
    # cluster sizes
    a = pm.constant(np.array([1., 1., 1.]))
    p = pm.Dirichlet('p', a=a, shape=k)
    # ensure all clusters have some points
    p_min_potential = pm.Potential('p_min_potential', tt.switch(tt.min(p) < .1, -np.inf, 0))


    # cluster centers
    means = pm.Normal('means', mu=[0, 1.5, 3], sd=1, shape=k)
    # break symmetry
    order_means_potential = pm.Potential('order_means_potential',
                                     tt.switch(means[1]-means[0] < 0, -np.inf, 0)
                                     + tt.switch(means[2]-means[1] < 0, -np.inf, 0))

    # measurement error
    sd = pm.Uniform('sd', lower=0, upper=2, shape=k)

    # …

我在使用 PYMC3 计算后验时遇到此错误：

with pm.Model() as model:
    p = pm.Gamma('p', alpha=1, beta=3, shape=regions.shape[0])
    q = pm.Gamma('q', alpha=1, beta=3, shape=regions.shape[0])
    m = pm.Lognormal('m', mu=np.log(total_M), sd=.25, shape=regions.shape[0])
    t = pm.Uniform('t', lower=0, upper=100, observed=sales.t)
    cid = pm.Categorical('cid', p=np.repeat(1./sales.shape[0], sales.shape[0]), observed=sales.region )
    sigma = pm.Gamma('sigma', alpha=1, beta=3)
    mu = m[cid]*(((p[cid]+q[cid])**2)/p[cid])*((np.exp(-(p[cid]+q[cid])*t))/((1+(q[cid]/p[cid])*np.exp(-(p[cid]+q[cid])*t))**2))
    Y_obs = pm.Normal('Ft', mu= mu, sd=sigma, observed= sales.sales)
    trace = pm.sample(100000,init = 'adapt_diag', progressbar = True, tune  = 1000) 

我尝试将 mu = mu 更改为 mu = np.log(mu)，它解决了错误，但与我的其他伙伴相比，结果却很糟糕。

我正在尝试运行以下模型，但在编译过程中失败：

import numpy as np
import pymc3 as pm


def sample_data(G=1, K=2):
    # mean proportion ([0,1]) for each g
    p_g = np.random.beta(2, 2, size=G)

    # concentration around each p_g
    c_g = np.random.lognormal(mean=0.5, sigma=1, size=G)

    # reparameterization for standard Beta(a,b)
    a_g = c_g * p_g / np.sqrt(p_g**2 + (1.-p_g)**2)
    b_g = c_g*(1.-p_g) / np.sqrt(p_g**2 + (1.-p_g)**2)

    # for each p_g, sample K proportions
    p_gk = np.random.beta(a_g[:, np.newaxis], b_g[:, np.newaxis], size=(G, K))

    return p_gk

# Data size
G = 3
K = …

我正在尝试将 PYMC3 用于贝叶斯模型，我想在新的看不见的数据上反复训练我的模型。我想我每次看到数据时都需要用先前训练过的模型的后验更新先验，类似于这里https://docs.pymc.io/notebooks/updating_priors.html 的实现方式。他们使用以下函数从样本中找到 KDE，并使用对 from_posterior 的调用替换模型中参数的每个原始定义。

def from_posterior(param, samples):
    smin, smax = np.min(samples), np.max(samples)
    width = smax - smin
    x = np.linspace(smin, smax, 100)
    y = stats.gaussian_kde(samples)(x)

    # what was never sampled should have a small probability but not 0,
    # so we'll extend the domain and use linear approximation of density on it
    x = np.concatenate([[x[0] - 3 * width], x, [x[-1] + 3 * width]])
    y = np.concatenate([[0], y, [0]])
    return Interpolated(param, x, y)

这是我的原始模型。 …

我是 pymc3 的新手，我读过《黑客贝叶斯方法》，并尽最大努力完成 pymc3 中现有的生存分析教程。但是，我不明白如何编写/解释“生存函数”。

对于这个问题，我从 NIST 定义的威布尔分布中生成了一些虚拟数据：

n = 1000
alpha = 1 
gam = 0.5
mu = 0
noise = np.random.normal(0, 0.025, [n, 1])

x = np.random.rand(n, 1)*10
f_x = (gam/alpha)*(((gam-mu)/alpha)**(gam-1))*np.exp(-((x-mu)/alpha)**gam)
y = f_x + noise

由于我想创建一个包含审查数据和未经审查数据的模型，例如贝叶斯参数化的 pymc3教程，因此我实现了截止并审查了这些数据点：

cens = np.array([1 if k < 7.5 else 0 for k in x])

然后我开始根据先验建立我的模型：

with pm.Model() as survival_model:

     alpha0 = pm.Normal('alpha0', mu=1, sigma = 1)
     gam0 = pm.Normal('gam0', mu=0.5, sigma = 1)
     mu0 = pm.Normal('mu0', mu=0.0, …

我试图在 Windows 中使用 Miniconda 安装 pymc3。安装当然不是一件简单的事，我花了很多时间进行故障排除才能了解发生了什么。

我读过很多不同人关于这个问题的帖子，发现本指南非常有帮助：

Windows 中 pymc3 安装指南

按照此处推荐的步骤进行操作后，我importing pymc3 as pm在 jupyter-notebook 中仍然遇到了警告：

WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.

以下是我已经测试过的一些解决方案，但无法消除此警告：

1. conda install -c conda forge libpython blas mkl-service m2w64-toolchain除了第一个 conda 环境中与 pymc3 (theano-pymc3和) 捆绑的软件包。arviz

2. 下载 GCC 编译器并确保它按照第一个指南和官方指南的建议以正确的顺序位于全局路径上在第二个 conda 环境中安装官方 pymc3m2w64-toolchain并且安装时没有.

3. 我还尝试使用 Rtools 中的 GCC 编译器，并更新了第三个 conda 环境中的环境变量中的路径。

我想我已经用尽了我读过的所有不同选项，但仍然无法删除此警告。在解决此警告问题时我是否遗漏了什么？

我知道这个警告不会停止计算，只会影响计算性能。因为我有相当大的数据集，所以我当然希望有一个可行的解决方案来解决这个问题:)

感谢您的任何反馈！