私は単純な PyMC 3 モデルを構築しようとしています。このモデルでは、潜在的な二変量ガウス密度で 2 つのカット ポイントと相関パラメーターを推定し、(多項) カウントのベクトルに対して 4 つの予測確率を生成します。(これは最終的に、これらのパラメータやその他のパラメータが多数の潜在的な多変量ガウス密度に対して推定される、より大きなモデルの一部になることを願っています。)
したがって、カット ポイント cx と cy を通常の確率変数としてモデル化し、相関パラメーター rho をスケーリングされたベータ確率変数としてモデル化したい (補足として、rho を処理するより良い方法を聞きたい - PyMC 3 はたとえば、通常の確率変数を切り捨てましたか?)。関数 mvnun を使用して、cx、cy、および rho の特定の値の予測確率を計算します。関数 mvnun は scipy.stats.mvn の一部です。これは、非常に正確な多変量正規 CDF 値を計算するための 2 つの関数を持つコンパイル済みの Fortran コードです。
相関行列 S に rho を貼り付けようとするか、cx または cy を積分限界を示す配列に入れると、次のようになります。
ValueError: setting an array element with a sequence.
cx、cy、および/または rho に固定数値を使用すると、mvnun は問題なく動作します ('with model:' ブロックの内外で)。PyMC RV がこのエラーを引き起こす理由を突き止めようとして、いろいろ調べてみましたが、困惑しています。cx、cy、および rho が theano TensorVariables であると私は収集しますが、theano TensorVariables について、これらの問題を引き起こす原因が何かを理解することはできません。
PyMC RV を引数として Fortran 関数を呼び出そうとすると、根本的な問題はありますか? または、私のコードに何らかの欠陥がありますか? 両方?まったく別の何か?
私はPyMCを初めて使用し、PyMC 3を現在のバージョンだと思ってインストールしました(数週間前にインストールしたときにアルファバージョンが存在しなかったことに注意してください)。おそらく、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)