1

Cam Davidson-Pilon のProbabilistic Programming & Bayesian Methods for Hackersを読んだ後、PyMC を使用して隠れマルコフ モデル (HMM) の学習問題に挑戦することにしました。今のところ、コードは連携していませんが、トラブルシューティングを通じて、問題の原因が絞り込まれたと感じています。

コードを小さなチャンクに分割し、t=0 での初期確率と放出確率に注目すると、時間 t=0 での単一状態の放出/観測パラメータを学習できます。しかし、別の状態 (合計 2 つの状態) を追加すると、データ入力に関係なく、パラメーター学習の結果は同じ (そして正しくない) になります。そのため、コードの一部で何か間違ったことをしたに違いないと感じています。これにより、初期確率関数@pm.deterministicからのサンプリングが許可されません。Init

コードのこの部分を使用して、状態 0 と状態 1 にそれぞれ対応する初期確率 p_bernと放出確率 p_0を学習することを目指しています。p_1放出は、私が自分の関数で表現しようとしている状態を条件としてい@pm.deterministicます。この決定的関数に「if」ステートメントを含めることはできますか? それが問題の根源のようです。

# This code is to test the ability to discern between two states with emissions

import numpy as np
import pymc as pm
from matplotlib import pyplot as plt

N = 1000
state = np.zeros(N)
data = np.zeros(shape=N)

# Generate data
for i in range(N):
    state[i] = pm.rbernoulli(p=0.3)
for i in range(N):
    if state[i]==0:
        data[i] = pm.rbernoulli(p=0.4)
    elif state[i]==1:
        data[i] = pm.rbernoulli(p=0.8)

# Prior on probabilities
p_bern = pm.Uniform("p_S", 0., 1.)
p_0 = pm.Uniform("p_0", 0., 1.)
p_1 = pm.Uniform("p_1", 0., 1.)

Init = pm.Bernoulli("Init", p=p_bern) # Bernoulli node

@pm.deterministic
def p_T(Init=Init, p_0=p_0, p_1=p_1, p_bern=p_bern):
    if Init==0:
        return p_0
    elif Init==1:
        return p_1

obs = pm.Bernoulli("obs", p=p_T, value=data, observed=True)
model = pm.Model([obs, p_bern, p_0, p_1])
mcmc = pm.MCMC(model)
mcmc.sample(20000, 10000)
pm.Matplot.plot(mcmc)

私はすでに次のことを無駄に試みました：

  1. @pm.potentialデコレーターを使用して共同ディストリビューションを作成する
  2. 自分の場所の配置を変更するInit(どこに配置すればよいかわからないコードのコメントを見ることができます)
  3. これ@pm.stochasticに似たものを使用してください

編集:クリスの提案に従って、ベルヌーイ ノードを決定論の外に移動しました。また、トラブルシューティングを容易にするために、コードをより単純なモデル (多項式ではなくベルヌーイ観測) に更新しました。

お時間をいただきありがとうございます。どんなフィードバックでも暖かく受け取られます。また、情報が不足している場合はお知らせください。

4

2 に答える 2

2

あなたが提供した更新された情報に基づいて、ここに動作するいくつかのコードがあります：

import numpy as np
import pymc as pm
from matplotlib import pyplot as plt

N = 1000
state = np.zeros(N)
data = np.zeros(shape=N)

# Generate data
state = pm.rbernoulli(p=0.3, size=N)
data = [int(pm.rbernoulli(0.8*s or 0.4)) for s in state]

# Prior on probabilities
p_S = pm.Uniform("p_S", 0., 1.)
p_0 = pm.Uniform("p_0", 0., 1.)
p_1 = pm.Uniform("p_1", 0., 1.)

# Use values of Init as indices to probabilities
Init = pm.Bernoulli("Init", p=p_S, size=N) # Bernoulli node

p_T = pm.Lambda('p_T', lambda p_0=p_0, p_1=p_1, i=Init: np.array([p_0, p_1])[i.astype(int)])

obs = pm.Bernoulli("obs", p=p_T, value=data, observed=True)
model = pm.MCMC(locals())
model.sample(20000, 10000)
model.summary()

データ生成ステップで、状態を使用して適切な真の確率にインデックスを付けていることに注意してください。の仕様でも基本的に同じことをしていp_Tます。これはかなりうまく機能しているように見えますが、どこで初期化されるかによって、 と の 2 つの値がいずれかの真の値に対応することに注意してくださいp_0(p_1一方が他方よりも大きくなるように制約するものは何もありません)。したがって、 の値はp_S真の状態確率の補数になる可能性があります。

于 2015-08-30T21:36:25.713 に答える