标签: pymc

我正在更新一些计算,其中我使用pymc2到pymc3,当我在模型上有一些离散的随机变量时,我在采样器行为方面遇到了一些问题.例如,请考虑使用pymc2的以下模型:

import pymc as pm

N = 100
data = 10

p = pm.Beta('p', alpha=1.0, beta=1.0)
q = pm.Beta('q', alpha=1.0, beta=1.0) 
A = pm.Binomial('A', N, p)
X = pm.Binomial('x', A, q, observed=True, value=data)

它并不代表任何东西,它只是一个模型,其中一个未观察到的变量是离散的.当我用pymc2对这个模型进行采样时,得到以下结果:

mcmc = pm.MCMC(model)
mcmc.sample(iter=100000, burn=50000, thin=100)
plot(mcmc)

但是当我尝试使用PYMC3时,我得到了这个:

with pm.Model() as model:
    N = 100
    p = pm.Beta('p', alpha=1.0, beta=1.0)
    q = pm.Beta('q', alpha=1.0, beta=1.0) 
    A = pm.Binomial('A', N, p)
    X = pm.Binomial('x', A, q, observed=10)

with model:
    start = pm.find_MAP()

with model:
    step = pm.NUTS() …

我试图结合cvxopt(优化求解器)和PyMC(采样器)来解决凸随机优化问题.

作为参考,安装两个软件包pip非常简单:

pip install cvxopt
pip install pymc 

两个包都可以很好地独立工作.以下是如何使用以下方法解决LP问题的示例cvxopt:

# Testing that cvxopt works
from cvxopt import matrix, solvers

# Example from http://cvxopt.org/userguide/coneprog.html#linear-programming

c = matrix([-4., -5.])
G = matrix([[2., 1., -1., 0.], [1., 2., 0., -1.]])
h = matrix([3., 3., 0., 0.])
sol = solvers.lp(c, G, h)
# The solution sol['x'] is correct: (1,1)

但是,当我尝试将它与PyMC一起使用时(例如,通过在其中一个系数上放置分布),PyMC会给出错误:

import pymc as pm
import cvxopt

c1 = pm.Normal('c1', mu=-4, tau=.5**-2)

@pm.deterministic
def my_lp_solver(c1=c1): …

我正在尝试在PyMC3中对多个链进行采样.在PyMC2中我会做这样的事情:

for i in range(N):
    model.sample(iter=iter, burn=burn, thin = thin)

我应该如何在PyMC3中做同样的事情？我看到'sample'方法中有一个'njobs'参数,但是当我为它设置一个值时它会抛出一个错误.我想使用采样链来获得'pymc.gelman_rubin'输出.

我之前实现了原始的贝叶斯概率矩阵分解(BPMF)模型pymc3.请参阅我之前的问题以供参考,数据源和问题设置.根据@twiecki对该问题的回答,我已经使用priors实现了模型的变体,LKJCorr用于相关矩阵和标准偏差的均匀先验.在原始模型中,协方差矩阵是从Wishart分布中提取的,但由于当前的局限性pymc3,Wishart分布无法从中正确采样.这个松散相关问题的答案为LKJCorr先验的选择提供了简洁的解释.新模型如下.

import pymc3 as pm
import numpy as np
import theano.tensor as t


n, m = train.shape
dim = 10  # dimensionality
beta_0 = 1  # scaling factor for lambdas; unclear on its use
alpha = 2  # fixed precision for likelihood function
std = .05  # how much noise to use for model initialization

# We will use separate priors for sigma and …

是否可以在pyMC3中逐步更新模型.我目前无法找到相关信息.所有文档始终使用先验已知数据.

但在我的理解中,贝叶斯模型也意味着能够更新信念.这在pyMC3中是否可行？我在哪里可以找到相关信息？

谢谢 :)

PyMC在其当前可用版本中如何适合用于连续发射HMM建模？

我对建立一个框架很感兴趣，在该框架中我可以轻松探索模型的变化，而不必更新E步骤和M步骤，并且无需为每次对模型所做的更改进行动态编程递归。

更具体的问题是：

• 在PyMC中为HMM建模时，我可以回答人们想解决的“典型”任务吗？即，除了参数估计之外，还可以推断出最可能的序列（通常用维特比算法完成），或者解决平滑问题？
• 与使用Expectation Maximization的实现相比，我希望基于采样的方法会更慢。如果这样可以使我在模型构建方面更具灵活性，那很好。我可以想象使用PyMC进行原型设计。不过，我想知道，是否可以期望PyMC在超过10k观测值的模型中在任何合理的时间内完成推断。
• 您是否建议从PyMC2或PyMC3开始进行模型构建。我知道推理引擎在版本之间会发生变化，因此我尤其想知道哪种类型的采样器可能更适合。

如果您认为PyMC在我的用例中不是一个很好的选择，那肯定也可以作为答案。

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

我发现使用 GPU 与 CPU 时 PYMC3 速度的差异为零。

我正在拟合一个需要 500K+ 样本才能收敛的模型。显然它非常慢，所以我尝试使用 GPU 来加快速度（在 EC2 上使用 GPU 实例）。Theano 报告正在使用 GPU，因此我相信 CUDA/Theano 配置正确。然而，我强烈怀疑 Pymc3 没有利用 GPU。

• 我是否需要显式将变量设置为 TensorType(float32, scalar) ？目前，它们是 float64。
• 是否只有部分采样器/可能性可以从 CUDA 中受益？我正在拟合基于泊松的模型，因此使用 Metropolis 采样器，而不是 NUTS
• 有没有办法检查 pymc3 是否正在使用 GPU？

我正在尝试将 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, …