标签: pymc

要学习PyMC,我正在尝试做一个简单的隐马尔可夫模型,如下所示:

with pymc3.Model() as hmm:
    # Transition "matrix"
    a_t = np.ones(num_states)
    T = [pymc3.Dirichlet('T{0}'.format(i), a = a_t,  shape = num_states) for i in xrange(num_states)]
    # Emission "matrix"
    a_e = np.ones(num_emissions)
    E = [pymc3.Dirichlet('E{0}'.format(i), a = a_e,  shape = num_emissions) for i in xrange(num_states)]
    # State models
    p0 = np.ones(num_states) / num_states
    # No shape, so each state is a scalar tensor
    states = [pymc3.Categorical('s0', p = p0)]
    emissions = [pymc3.Categorical('z0',
                               p = ifelse(eq(states[0], 0), E[0], ifelse(eq(states[0], 1), E[1], E[2])),
                               observed …

数据集中有两列,分别是user_id和site_name.它记录每个用户浏览的每个站点名称.

toy_dict = {'site_name': {0: u'\u4eac\u4e1c\u7f51\u4e0a\u5546\u57ce',
1: u'\u963f\u91cc\u4e91',
2: u'\u6dd8\u5b9d\u7f51',
3: u'\u624b\u673a\u6dd8\u5b9d\u7f51',
4: u'\u6211\u4eec\u7684\u70b9\u5fc3\u7f51',
5: u'\u8c46\u74e3\u7f51',
6: u'\u9ad8\u5fb7\u5730\u56fe',
7: u'\u817e\u8baf\u7f51',
8: u'\u70b9\u5fc3',
9: u'\u767e\u5ea6',
10: u'\u641c\u72d7',
11: u'\u8c37\u6b4c',
12: u'AccuWeather\u6c14\u8c61\u9884\u62a5',
13: u'\u79fb\u52a8\u68a6\u7f51',
14: u'\u817e\u8baf\u7f51',
15: u'\u641c\u72d7\u7f51',
16: u'360\u624b\u673a\u52a9\u624b',
17: u'\u641c\u72d0',
18: u'\u767e\u5ea6'},
'user_id': {0: 37924550,
1: 37924550,
2: 37924550,
3: 37924550,
4: 37924550,
5: 37924550,
6: 37924550,
7: 37924550,
8: 37924551,
9: 37924551,
10: 37924551,
11: 37924551,
12: 37924551,
13: 37924552,
14: 45285152,
15: 45285153,
16: 45285153,
17: 45285153,
18: …

假设我们在X上给出了先验(例如X~高斯)和前向算子y = f(x).假设我们通过实验进一步观察到y,并且该实验可以无限期地重复.假设输出Y是高斯(Y~高斯)或无噪声(Y~Delta(观察)).

如何根据观察结果不断更新我们对X的主观知识水平？我用PyMC尝试了以下模型,但似乎我遗漏了一些东西:

from pymc import *

xtrue = 2                        # this value is unknown in the real application
x = rnormal(0, 0.01, size=10000) # initial guess

for i in range(5):
    X = Normal('X', x.mean(), 1./x.var())
    Y = X*X                        # f(x) = x*x
    OBS = Normal('OBS', Y, 0.1, value=xtrue*xtrue+rnormal(0,1), observed=True)
    model = Model([X,Y,OBS])
    mcmc = MCMC(model)
    mcmc.sample(10000)

    x = mcmc.trace('X')[:]       # posterior samples

后部没有收敛到xtrue.

我在这里发布了一个python笔记本:http://nbviewer.ipython.org/gist/awellis/9067358

我正在尝试使用PyMC 3创建一个probit回归模型,使用生成的数据来恢复已知参数(参见笔记本).拦截的估计几乎没有,但是斜率估计是偏离标记的.

我的模型看起来像这样:

with pm.Model() as model:

    # priors
    alpha = pm.Normal('alpha', mu=0, tau=0.001)
    beta = pm.Normal('beta', mu=0, tau=0.001)

    # linear predictor
    theta_p = (alpha + beta * x)

    # logic transform (just for comparison - this seems to work ok)
#     def invlogit(x):
#         import theano.tensor as t
#         return t.exp(x) / (1 + t.exp(x))
#     theta = invlogit(theta_p)


    # Probit transform: this doesn't work
    def phi(x):
        import theano.tensor as t
        return 0.5 * (1 …

当尝试通过conda安装pymc时,我收到以下内容:

C:\ Anaconda> conda install -c https://conda.binstar.org/pymc pymc获取包元数据:...错误:找不到匹配的包:pymc

安装来自pymc分发页面:https://binstar.org/pymc/pymc

我目前的anaconda版本是最新的:

C:\ Anaconda> conda update --prefix C:\ Anaconda anaconda获取包元数据:..解决包规范:.#已安装所有请求的包.

在C:\ Anaconda环境中的#包:

#

anaconda 1.9.2 np18py27_0

所以作为康达新手,我不太清楚我错过了什么.也许我必须首先授权binstar？(我相信没有代理问题.)

非常感谢您的建议!

我是PyMC的新手,并尝试使用最大后验估计来拟合我的非齐次泊松过程和分段恒定速率函数.

我的过程描述了一天中的一些事件.因此,我将一天分成24小时,这意味着,我的速率函数中有24个常数(分段常数).

结合以下思路:

• 黑客的概率规划和贝叶斯方法:第1章
• SO问题:使用可变速率拟合带顶点的泊松过程

我提出了以下一段代码,这是不满意的(结果明智,我确定这是错的):

import numpy as np
import pymc

eventCounter = np.zeros(24) # will be filled with real counts before going on
alpha = 1.0 / eventCounter.mean()

a0 = pymc.Exponential('a0', alpha)
a1 = pymc.Exponential('a1', alpha)
a2 = pymc.Exponential('a2', alpha)
a3 = pymc.Exponential('a3', alpha)
a4 = pymc.Exponential('a4', alpha)
a5 = pymc.Exponential('a5', alpha)
a6 = pymc.Exponential('a6', alpha)
a7 = pymc.Exponential('a7', alpha)
a8 = pymc.Exponential('a8', alpha)
a9 = pymc.Exponential('a9', alpha)
a10 = pymc.Exponential('a10', alpha) …

如何在PyMC3中定义自定义可能性？在PyMC2中,我可以使用@pymc.potential.我试图用pymc.Potential在PyMC3,但是,似乎布尔操作不能被应用到的参数(我得到这样一个错误这样,当我做到这一点).例如,以下代码不起作用:

from pymc import *

with Model() as model:
    x = Normal('x', 1, 1)

    def z(u):
        if u > 0: #comparisons like this are not supported
        # if theano.tensor.lt(0,u): this is how comparison should be done
            return u ** 2
        return -u**3

    x2 = Potential('x2', z(x))

    start = model.test_point
    h = find_hessian(start)
    step = Metropolis(model.vars, h)
    sample(100, step, start)

我不可能将可能性内的所有比较更改为Theano语法(即theano.tensor.{lt,le,eq,neq,gt,ge}).无论如何使用定义类似于PyMC2的似然函数？

我仍在学习PYMC3，但是我在文档中找不到以下问题的任何内容。考虑这个问题的贝叶斯结构时间序列（BSTS）模型，没有季节性。可以在PYMC3中对其进行建模，如下所示：

import pymc3, numpy, matplotlib.pyplot
# generate some test data
t = numpy.linspace(0,2*numpy.pi,100)
y_full = numpy.cos(5*t)
y_train = y_full[:90]
y_test = y_full[90:]

# specify the model
with pymc3.Model() as model:
  grw = pymc3.GaussianRandomWalk('grw',mu=0,sd=1,shape=y_train.size)
  y = pymc3.Normal('y',mu=grw,sd=1,observed=y_train)
  trace = pymc3.sample(1000)
  y_mean_pred = pymc3.sample_ppc(trace,samples=1000,model=model)['y'].mean(axis=0)

  fig = matplotlib.pyplot.figure(dpi=100)
  ax = fig.add_subplot(111)
  ax.plot(t,y_full,c='b')
  ax.plot(t[:90],y_mean_pred,c='r')
  matplotlib.pyplot.show()

现在，我想预测接下来10个时间步的行为，即y_test。我还想在此区域包括产生贝叶斯锥的可靠区域，例如，请参见此处。不幸的是，在上述连接中产生锥体的机构有点模糊。在更传统的AR模型中，人们可以学习平均回归系数并手动扩展平均曲线。但是，在此BSTS模型中，没有明显的方法来执行此操作。或者，如果存在回归变量，则可以使用theano.shared并使用更细/扩展的网格对其进行更新，以使用sample_ppc进行插补和推断，但这在此设置中并不是真正的选择。也许sample_ppc在这里是个红鲱鱼，但尚不清楚如何进行其他操作。任何帮助都将受到欢迎。

我试图将简单的生存模型从这里(介绍中的第一个)从PyMC 2移植到PyMC 3.然而,我没有找到任何等同于"观察"的装饰器,并且我尝试编写新的分发失败了.有人可以提供一个例子,如何在PyMC 3中完成这项工作？

我是比较新的PYMC3，我试图实现贝叶斯结构时间序列（BSTS）没有回归系数，比如模型拟合这里在R.模型如下：

我可以使用GaussianRandomWalk来实现局部线性趋势，如下所示：

delta = pymc3.GaussianRandomWalk('delta',mu=0,sd=1,shape=99)
mu = pymc3.GaussianRandomWalk('mu',mu=delta,sd=1,shape=100)

但是，我对如何在PYMC3中编码季节性变量（tau）感到困惑。我是否需要参加自定义的随机游走课程，还是有其他技巧？