我正在尝试构建一个简单的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：'块中）或者在'with model：'块之内 . 我一直在四处寻找，试图找出PyMC RV导致这个错误的原因，但我很难过 . 我认为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)
1 回答
这在PyMC 2.3中会稍微容易一些，因此安装2.3是一个不错的选择 .
在PyMC 3中不起作用，因为cx，cy和rho是theano TensorVariables而mvncdf需要numpy数组 . Theano变量不是一种numpy数组 . 相反，它们代表了导致numpy数组的计算 .
你现在是一个非常简单的方法，尽管it is coming . 您还可以查看guide to making a Theano Op .