0

PyMC を使用してロジスティック回帰モデルを解こうとしました。ただし、診断プロットは非常に高い自己相関を示し、事後分布から繰り返しサンプリングした後、非常に異なる結果が得られることがあるため、PyMC を適切に使用していない可能性があります。

モデルは次のとおりです。

Y_i = Bernoulli(p_i)
logit(p_i) = b0 + b1*x1 + b2*x2

実装は次のとおりです。

import pymc as pm
import numpy as np

b0 =  pm.Normal('b0',  0., 1e-6, value=0.)
b1 =  pm.Normal('b1',  0., 1e-6, value=0.)
b2 =  pm.Normal('b2',  0., 1e-6, value=0.)

x1 = np.array([31, 31, 36, 30, 32, 33, 31, 33, 32])
x2 = np.array([31, 16, 35, 42, 19, 37, 29, 23, 15])
value = np.array([1, 1, 1, 1, 0, 1, 1, 1, 0])

@pm.deterministic
def p(b0=b0, b1=b1, b2=b2, x1=x1, x2=x2):
    return np.exp(b0 + b1*x1 + b2*x2)/(1.0 + np.exp(b0 + b1*x1 + b2*x2))

obs = pm.Bernoulli('obs', p=p, value=value, observed=True)

m = pm.MCMC([obs, b0, b1, b2])

サンプリングしm.sample(500000, 200000, 50)て結果の事後分布をプロットすると、次のようになります。

ここに画像の説明を入力 ここに画像の説明を入力 ここに画像の説明を入力

より良い結果を得るための 2 回目の試みでは、@pm.observed を使用しました。

import pymc as pm
import numpy as np

b0 =  pm.Normal('b0',  0., 1e-6, value=0.)
b1 =  pm.Normal('b1',  0., 1e-6, value=0.)
b2 =  pm.Normal('b2',  0., 1e-6, value=0.)

x1 = np.array([31, 31, 36, 30, 32, 33, 31, 33, 32])
x2 = np.array([31, 16, 35, 42, 19, 37, 29, 23, 15])
value = np.array([1, 1, 1, 1, 0, 1, 1, 1, 0])

p = b0 + b1*x1 + b2*x2

@pm.observed
def obs(p=p, value=value):
    return pm.bernoulli_like(value, pm.invlogit(p))

m = pm.MCMC([obs, p, b0, b1, b2])

しかし、それは高い自己相関も生み出します。

サンプルサイズを増やしましたが、あまり成功しませんでした。私は何が欠けていますか？

4

1 に答える 1

2

あなたは間違いなくいくつかの収束の問題を抱えています。問題の一部は、サンプル サイズが非常に小さい (n=9) にもかかわらず、3 つのパラメーターを適合させようとしていることです。パラメータをブロック更新することでマイレージは少し良くなりますが、インターセプトはまだうまくいきません:

beta =  pm.Normal('beta',  0., 0.001, value=[0.]*3)

x1 = np.array([31, 31, 36, 30, 32, 33, 31, 33, 32])
x2 = np.array([31, 16, 35, 42, 19, 37, 29, 23, 15])
value = np.array([1, 1, 1, 1, 0, 1, 1, 1, 0])

p = pm.Lambda('p', lambda b=beta: pm.invlogit(b[0] + b[1]*x1 + b[2]*x2))

obs = pm.Bernoulli('obs', p=p, value=value, observed=True)

(ところで、パラメーターをロジット変換する場合、ベータ版で超拡散事前確率は必要ありません。0.001 より小さいとかなり無駄になります)。ベータ版でアダプティブ Metropolis サンプラーを使用することで、さらに改善されます。

M.use_step_method(AdaptiveMetropolis, M.beta)

これにより収束しますが、ご覧のとおり、事前に通知する情報はありません。

まとめ

于 2014-09-01T20:09:43.033 に答える