私は PyMC を使用して、病状の 4 つの異なる予測因子を入力として取り、それらを組み合わせて、「はい、この患者にはこの状態があります」という予測因子のサブセットが与えられた場合に、患者がその状態にあるという全体的な事後確率を生成しようとしています。
アイデアは、ベータ分布から各予測子のシータ (条件の全体的な割合) と偽陰性および偽陽性の割合を選択し、そこからベイズの定理を使用して周辺確率と事後確率を計算することです。カウントの 16 個の観測値の配列があり、予測子の可能な組み合わせごとに 1 つです (4 つの予測子があるため、予測子の可能な組み合わせは 2**4 = 16 通りあります)。PyMC チュートリアルhttp://pymc-devs.github.io/pymc/tutorialの次の例で、disasters_array をポアソン分布で使用する方法と同様に、この最後のカウント セットと限界確率を多項分布に入力しています。 .html .
これを試して実行するために私が書いたコードは次のとおりです。
from pymc import Multinomial, Beta, deterministic
from numpy import array
n = 4 # number of predictors
counts_array = array([2942808473, 17491655, 21576, 23189, 339805, 89159, 168214, 76044, 43138288, 530963, 22682, 22169, 462052, 129472, 2804257, 3454104]) # 16 counts - one count value for each possible permutation of predictors that detected medical condition
pred = array([[0,0,0,0],[1,0,0,0],[0,1,0,0],[1,1,0,0],[0,0,1,0],[1,0,1,0],[0,1,1,0],[1,1,1,0],[0,0,0,1],[1,0,0,1],[0,1,0,1],[1,1,0,1],[0,0,1,1],[1,0,1,1],[0,1,1,1],[1,1,1,1]]); # array of yes/no's from predictors for each value in counts_array above
theta = Beta('theta',1,2)
fn = fp = tn = tp = [0] * 4;
for i in range(0,n):
fn[i] = Beta('fn' + str(i),1,2)
fp[i] = Beta('fp' + str(i),1,2)
tn[i] = 1 - fp[i]
tp[i] = 1 - fn[i]
@deterministic(plot=False)
def margprobs():
mp = [0] * 2**n; # initialize with vector of 2**n zeros
for i in range(0,2**n):
mp[i] = theta *\
(fn[0]**(1-pred[i][0])) * (tp[0]**pred[i][0]) *\
(fn[1]**(1-pred[i][1])) * (tp[1]**pred[i][1]) *\
(fn[2]**(1-pred[i][2])) * (tp[2]**pred[i][2]) *\
(fn[3]**(1-pred[i][3])) * (tp[3]**pred[i][3])\
+ (1-theta) *\
(tn[0]**(1-pred[i][0])) * (fp[0]**pred[i][0]) *\
(tn[1]**(1-pred[i][1])) * (fp[1]**pred[i][1]) *\
(tn[2]**(1-pred[i][2])) * (fp[2]**pred[i][2]) *\
(tn[3]**(1-pred[i][3])) * (fp[3]**pred[i][3]);
return mp;
@deterministic(plot=False)
def postprobs():
pp = [0] * 2**n; # initialize with vector of 2**n zeros
for i in range(0,2**n):
pp[i] = theta *\
(fn[0]**(1-pred[i][0])) * (tp[0]**pred[i][0]) *\
(fn[1]**(1-pred[i][1])) * (tp[1]**pred[i][1]) *\
(fn[2]**(1-pred[i][2])) * (tp[2]**pred[i][2]) *\
(fn[3]**(1-pred[i][3])) * (tp[3]**pred[i][3])\
/ margprobs[i];
return pp;
counts = Multinomial(name="counts", value=counts_array, n=2**n, p=margprobs, observed=True)
これを実行すると、カウントを計算するときに、最後の行でエラーが発生します。
$ python test.py
Traceback (most recent call last):
File "test.py", line 46, in <module>
counts = Multinomial(name="counts", value=counts_array, n=2**n, p=margprobs, observed=True)
File "/Users/jtr4v/anaconda/lib/python2.7/site-packages/pymc/distributions.py", line 3269, in __init__
verbose=verbose, **kwds)
File "/Users/jtr4v/anaconda/lib/python2.7/site-packages/pymc/PyMCObjects.py", line 772, in __init__
if not isinstance(self.logp, float):
File "/Users/jtr4v/anaconda/lib/python2.7/site-packages/pymc/PyMCObjects.py", line 929, in get_logp
raise ZeroProbability(self.errmsg)
pymc.Node.ZeroProbability: Stochastic counts's value is outside its support,
or it forbids its parents' current values
明らかに、このエラーは、PyMC が Multinomial() に入力している値を気に入らないことに関係していますが、どちらが間違っているのかわかりません。値は counts_array (カウントの観測値) である必要があり、n は 16 である必要があります。これは、カウントに 16 項目の配列を選択し、予測子の可能な組み合わせごとに 1 つ、p が限界確率であり、観測されたものである必要があるためです。私は値を観察したので、真でなければなりません。
私は何を間違っていますか？
編集:それが役立つ場合、以前は次のコードを使用してR2jagsでこれを行っていました:
model {
theta ~ dbeta(1,2); # draw theta from a beta distribution
for (i in 1:N) { # draw an false positive and false negative rate for each of N predictors
fp[i] ~ dbeta(1,2);
fn[i] ~ dbeta(1,2);
tp[i] <- 1-fn[i]; # true positive = 1 - false negative rate
tn[i] <- 1-fp[i]; # true negative rate = 1 - false positive rate
}
for (j in 1:M) {
# Bayes theorem, i.e.
# posterior probability =
# P(A) * P(B|A) /
# /
# P(A) * P(B|A) + P(-A) * P(B|-A) # <--- last line is marginal probability
#
# x is a vector of 1's and 0's indicating whether the ith predictor said yes or no
margprobs[j] <- (theta *
(fn[1]^(1-x[j,1])) * (tp[1]^x[j,1]) *
(fn[2]^(1-x[j,2])) * (tp[2]^x[j,2]) *
(fn[3]^(1-x[j,3])) * (tp[3]^x[j,3]) *
(fn[4]^(1-x[j,4])) * (tp[4]^x[j,4])
+ (1-theta) *
(tn[1]^(1-x[j,1])) * (fp[1]^x[j,1]) *
(tn[2]^(1-x[j,2])) * (fp[2]^x[j,2]) *
(tn[3]^(1-x[j,3])) * (fp[3]^x[j,3]) *
(tn[4]^(1-x[j,4])) * (fp[4]^x[j,4]));
postprobs[j] <- theta *
(fn[1]^(1-x[j,1])) * (tp[1]^x[j,1]) *
(fn[2]^(1-x[j,2])) * (tp[2]^x[j,2]) *
(fn[3]^(1-x[j,3])) * (tp[3]^x[j,3]) *
(fn[4]^(1-x[j,4])) * (tp[4]^x[j,4])
/ margprobs[j];
}
counts ~ dmulti(margprobs, total);
}