PyMC 3 中的确定性变量和 Fortran Scipy 函数答案

【问题标题】：Deterministic variables and a Fortran Scipy function in PyMC 3PyMC 3 中的确定性变量和 Fortran Scipy 函数
【发布时间】：2014-04-02 00:22:38
【问题描述】：

我正在尝试构建一个简单的 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)

【问题讨论】：

标签： pymc theano pymc3

【解决方案1】：

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

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

您需要将 mvncdf 变成 Theano 运算符（numpy 函数的 theano 模拟）。不幸的是，尽管it is coming，目前还没有一种超级简单的方法可以做到这一点。你也可以看看guide to making a Theano Op。

【讨论】：