标签: pymc3

我遇到了一个常见的问题,我想知道是否有人可以提供帮助.我经常想在两种模式下使用pymc3:训练(即实际运行对参数的推断)和评估(即使用推断参数来生成预测).

总的来说,我想要一个后验过度的预测,而不仅仅是逐点估计(这是贝叶斯框架的一部分,不是吗？).当您的训练数据得到修复时,通常可以通过将类似形式的模拟变量添加到观察变量来实现.例如,

from pymc3 import *

with basic_model:

    # Priors for unknown model parameters
    alpha = Normal('alpha', mu=0, sd=10)
    beta = Normal('beta', mu=0, sd=10, shape=2)
    sigma = HalfNormal('sigma', sd=1)

    # Expected value of outcome
    mu = alpha + beta[0]*X1 + beta[1]*X2

    # Likelihood (sampling distribution) of observations
    Y_obs = Normal('Y_obs', mu=mu, sd=sigma, observed=Y)
    Y_sim = Normal('Y_sim', mu=mu, sd=sigma, shape=len(X1))

    start = find_MAP()
    step = NUTS(scaling=start)
    trace = sample(2000, step, start=start)

但是如果我的数据发生了变化呢 假设我想根据新数据生成预测,但不再重复推断.理想情况下,我有一个类似的功能predict_posterior(X1_new, X2_new, 'Y_sim', trace=trace),甚至predict_point(X1_new, X2_new, 'Y_sim', …

我一直在试图实现贝叶斯线性回归使用模型PyMC3与真实数据(即不是从线性函数+高斯噪声)从数据集sklearn.datasets.我选择了具有最小数量的属性(即load_diabetes())形状为的回归数据集(442, 10); 就是,442 samples和10 attributes.

我相信我的模型正在运行,后面看起来还不错,可以预测并弄清楚这些东西是如何起作用的......但我意识到我不知道如何使用这些贝叶斯模型进行预测!我试图避免使用glm和patsy符号,因为我很难理解使用它时实际发生了什么.

我尝试了以下内容: 从pymc3 和http://pymc-devs.github.io/pymc3/posterior_predictive/中的推断参数生成预测,但我的模型在预测时非常糟糕,或者我做错了.

如果我实际上正在做正确的预测(我可能不是),那么任何人都可以帮助我优化我的模型.我不知道是否最少mean squared error,absolute error或类似的东西在贝叶斯框架中有效.理想情况下,我想得到一个number_of_rows数组=我的X_te属性/数据测试集中的行数,以及来自后验分布的样本列数.

import pymc3 as pm
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
from scipy import stats, optimize
from sklearn.datasets import load_diabetes
from sklearn.cross_validation import train_test_split
from theano …

我刚刚完成了Osvaldo Martin的Python书中的贝叶斯分析(理解贝叶斯概念和一些花哨的numpy索引的好书).

我真的想将我的理解扩展到贝叶斯混合模型,用于无监督的样本聚类.我所有的谷歌搜索都让我看到了Austin Rochford的教程,这本教程非常有用.我理解发生了什么,但我不清楚它如何适应群集(特别是使用群集分配的多个属性,但这是一个不同的主题).

我知道如何分配先验,Dirichlet distribution但我无法弄清楚如何获得集群PyMC3.看起来大多数mus会聚到质心(即我从中采样的分布方式),但它们仍然是分开的components.我考虑过weights(w在模型中)截止,但这似乎不像我想象的那样工作,因为多个components具有稍微不同的平均参数mus正在收敛.

如何从此PyMC3模型中提取聚类(质心)？我给了它最多的15组件,我想收敛3.在mus似乎是在正确的位置,但权重搞砸b他们被其他集群之间分配/ C,所以我不能用一个权重阈值(除非我把它们合并,但我不认为这是事情是这样的通常做完).

import pymc3 as pm
import numpy as np
import matplotlib.pyplot as plt
import multiprocessing
import seaborn as sns
import pandas as pd
import theano.tensor as tt
%matplotlib inline

# Clip at 15 components
K = 15 …

有没有办法加快这个简单的PyMC模型？在20-40个数据点上,需要约5-11秒才能适应.

import pymc
import time
import numpy as np
from collections import OrderedDict

# prior probability of rain
p_rain = 0.5
variables = OrderedDict()
# rain observations
data = [True, True, True, True, True,
        False, False, False, False, False]*4
num_steps = len(data)
p_rain_given_rain = 0.9
p_rain_given_norain = 0.2
p_umbrella_given_rain = 0.8
p_umbrella_given_norain = 0.3
for n in range(num_steps):
    if n == 0:
        # Rain node at time t = 0
        rain = pymc.Bernoulli("rain_%d" %(n), p_rain)
    else:
        rain_trans = \
          pymc.Lambda("rain_trans", …

据我了解，Python 中使用上下文管理器来定义对象的初始化和完成代码段 (__enter__和)。__exit__

然而，在PyMC3 教程中，他们展示了以下上下文管理器示例：

basic_model = pm.Model()

with basic_model:

    # Priors for unknown model parameters
    alpha = pm.Normal('alpha', mu=0, sd=10)
    beta = pm.Normal('beta', mu=0, sd=10, shape=2)
    sigma = pm.HalfNormal('sigma', sd=1)

    # Expected value of outcome
    mu = alpha + beta[0]*X1 + beta[1]*X2

    # Likelihood (sampling distribution) of observations
    Y_obs = pm.Normal('Y_obs', mu=mu, sd=sigma, observed=Y)

并提到这的目的是将变量alpha、beta、sigma和mu与Y_obs模型 basic_model 相关联。

我想了解这样的机制是如何运作的。在我发现的上下文管理器的解释 中，我没有看到任何建议上下文块中定义的变量或对象如何以某种方式“关联”到上下文管理器。库（PyMC3）似乎可以以某种方式访问​​“当前”上下文管理器，因此它可以在幕后将每个新创建的语句与其关联。但是库如何访问上下文管理器呢？

有人可以提供一些关于如何并行化PyMC MCMC代码的一般性说明.我试图LASSO按照这里给出的例子运行回归.我在某处读到默认情况下并行采样,但是我是否还需要使用类似的功能Parallel Python来使其工作？

这是一些我希望能够在我的机器上并行化的参考代码.

x1 = norm.rvs(0, 1, size=n)
x2 = -x1 + norm.rvs(0, 10**-3, size=n)
x3 = norm.rvs(0, 1, size=n)

X = np.column_stack([x1, x2, x3])
y = 10 * x1 + 10 * x2 + 0.1 * x3

beta1_lasso = pymc.Laplace('beta1', mu=0, tau=1.0 / b)
beta2_lasso = pymc.Laplace('beta2', mu=0, tau=1.0 / b)
beta3_lasso = pymc.Laplace('beta3', mu=0, tau=1.0 / b)

@pymc.deterministic
def y_hat_lasso(beta1=beta1_lasso, beta2=beta2_lasso, beta3=beta3_lasso, x1=x1, x2=x2, x3=x3):
    return beta1 * x1 …

一些简单的theano代码完美运行,当我导入pymc3时停止工作

这里有一些snipets重现错误:

#Initial Theano Code (this works)    
import theano.tensor as tsr

x = tsr.dscalar('x')
y = tsr.dscalar('y')
z = x + y

#Snippet 1
import pymc3 as pm
import theano.tensor as tsr

x = tsr.dscalar('x')
y = tsr.dscalar('y')
z = x + y

#Snippet 2
import theano.tensor as tsr
import pymc3 as pm

x = tsr.dscalar('x')
y = tsr.dscalar('y')
z = x + y

#Snippet 3
import pymc3 as pm

x = pm.theano.tensor.dscalar('x')
y = pm.theano.tensor.dscalar('y')
z = x + …

在pymc3中如何配置截断的正常先验？在pymc2中,它非常简单(下图),但在pymc3中似乎不再有截断的正态分布.

Pymc2:

TruncatedNormal('gamma_own_%i_' % i, mu=go, tau=v_gamma_inv, value=0, a=-np.inf, b=0)

Pymc3 :?

我正在学习PyMC3进行贝叶斯建模.您可以使用以下命令创建模型和示例:

import pandas as pd
import pymc3 as pm

# obs is a DataFrame with a single column, containing
# the observed values for variable height
obs = pd.DataFrame(...)

# we create a pymc3 model
with pm.Model() as m:
    mu = pm.Normal('mu', mu=178, sd=20)
    sigma = pm.Uniform('sigma', lower=0, upper=50)
    height = pm.Normal('height', mu=mu, sd=sigma, observed=obs)
    trace = pm.sample(1000, tune=1000)

pm.traceplot(trace)

当我检查trace(在这种情况下来自后验概率的1000个样本)时,我注意到创建了2个链:

>>> trace.nchains
2

我阅读了关于PyMC3的教程,并查看了API,但我不清楚链表示什么(在这种情况下,我从后面询问了1000个样本,但我得到了2个链,每个链有1000个来自后面的样本).

链条是否具有相同参数的采样器的不同运行,或者它们是否具有其他一些含义/目的？

我正在努力在 pymc3 中实现具有自定义可能性的线性回归。

我之前在 CrossValidated 上发布过这个问题，建议在此处发布，因为该问题更面向代码（此处已关闭帖子）

假设您有两个自变量 x1、x2 和一个目标变量 y，以及一个称为 delta 的指示变量。

• 当 delta 为 0 时，似然函数为标准最小二乘
• 当 delta 为 1 时，仅当目标变量大于预测时，似然函数才是最小二乘贡献

观察到的数据的示例片段：

x_1  x_2     observed_target  
10    1   0   100              
20    2   0   50               
5    -1   1   200             
10   -2   1   100             

有谁知道如何在 pymc3 中实现这一点？作为起点...

model =  pm.Model()
with model as ttf_model:

  intercept = pm.Normal('param_intercept', mu=0, sd=5)
  beta_0 = pm.Normal('param_x1', mu=0, sd=5)
  beta_1 = pm.Normal('param_x2', mu=0, sd=5)
  std = pm.HalfNormal('param_std', beta = 0.5)

  x_1 = pm.Data('var_x1', df['x1'])
  x_2 …