1

我正在尝试构建一个简单的 PyMC 3 模型，在该模型中，我估计潜在二元高斯密度中的两个切点和一个相关参数，为（多项式）计数向量产生四个预测概率。（我希望这最终会成为一个更大模型的一部分，在这个模型中，这些参数和其他参数是针对许多潜在的多元高斯密度进行估计的。）

因此，我想将切点 cx 和 cy 建模为正常随机变量，并将相关参数 rho 建模为缩放的 Beta 随机变量（作为旁注，我很想听到处理 rho 的更好方法 - PyMC 3截断了正常的随机变量，例如？）。我想使用函数 mvnun 来计算给定 cx、cy 和 rho 值的预测概率。函数 mvnun 是 scipy.stats.mvn 的一部分，它是一段编译的 Fortran 代码，有两个函数用于计算非常准确的多元正态 CDF 值。

如果我尝试将 rho 粘贴在相关矩阵 S 中，或者如果我将 cx 或 cy 放入指示积分限制的数组中，我会得到：

ValueError: setting an array element with a sequence.

如果我对 cx、cy 和/或 rho 使用固定数值，mvnun 就可以正常工作（在“with model:”块内或外）。我一直在四处寻找，试图弄清楚为什么 PyMC RVs 会导致这个错误，但我很难过。我收集到 cx、cy 和 rho 是 theano TensorVariables，但我无法弄清楚关于 theano TensorVariables 的什么（如果有的话）会导致这些问题。

尝试使用 PyMC RV 作为参数调用 Fortran 函数是否存在根本问题？还是我的代码在某些方面存在缺陷？两个都？完全不同的东西？

我是 PyMC 的新手，我安装 PyMC 3 时认为它是当前版本（我发誓说几周前我安装它时没有 alpha 版本）。也许我应该安装 2.3 并弄清楚如何把它和那个放在一起？

无论如何，任何有关如何解决问题的建议都将不胜感激。

这是我的代码：

import pymc as pm
import numpy as np
from scipy.stats.mvn import mvnun as mvncdf

counts = np.array([100,35,45,20])
N = np.sum(counts)

def scaleBeta(x):
    return x*1.98 - 0.99

model = pm.Model()

with model:

    cx = pm.Normal('Cx', mu=0., tau=1.)
    cy = pm.Normal('Cy', mu=0., tau=1.)

    mu = np.array([0., 0.])

    rho_beta = pm.Beta('rho_beta', alpha=2, beta=2)
    rho = pm.Deterministic('rho',scaleBeta(rho_beta))
    S = np.array([1, rho, rho, 1]).reshape(2,2)

    low_aa = np.array([-50.,-50.])
    upp_aa = np.array([cx, cy])
    low_ab = np.array([-50., cy])
    upp_ab = np.array([cx, 50.])
    low_ba = np.array([cx, -50.])
    upp_ba = np.array([50., cy])
    low_bb = np.array([cx, cy])
    upp_bb = np.array([50., 50.])

    p_aa,i = mvncdf(lower=low_aa,upper=upp_aa,means=mu,covar=S)
    p_ab,i = mvncdf(lower=low_ab,upper=upp_ab,means=mu,covar=S)
    p_ba,i = mvncdf(lower=low_ba,upper=upp_ba,means=mu,covar=S)
    p_bb,i = mvncdf(lower=low_bb,upper=upp_bb,means=mu,covar=S)

    pv = np.array([p_aa,p_ab,p_ba,p_bb])

    cv = pm.Multinomial('counts',n=N,p=pv,observed=counts)
4

1 回答 1

1

这在 PyMC 2.3 中会稍微容易一些，因此安装 2.3 是一个不错的选择。

在 PyMC 3 中不起作用，因为 cx、cy 和 rho 是 theano TensorVariables，而 mvncdf 需要 numpy 数组。Theano 变量不是一种 numpy 数组。相反，它们表示导致 numpy 数组的计算。

您需要将 mvncdf 变成 Theano 运算符（theano 模拟 numpy 函数）。不幸的是，目前还没有一种超级简单的方法可以做到这一点，尽管它即将到来。您还可以查看制作 Theano Op 的指南。

于 2014-04-02T06:54:31.207 回答