1

これはPyMC のフォローアップです: マルコフ システムでのパラメーター推定

各タイムステップでの位置と速度によって定義されるシステムがあります。システムの動作は次のように定義されます。

vel = vel + damping * dt
pos =  pos + vel * dt

これが私の PyMC モデルです。見積もりvel、posそして最も重要なことdamping。

# PRIORS
damping = pm.Normal("damping", mu=-4, tau=(1 / .5**2))
# we assume some system noise
tau_system_noise = (1 / 0.1**2)

# the state consist of (pos, vel); save in lists
# vel: we can't judge the initial velocity --> assume it's 0 with big std
vel_states = [pm.Normal("v0", mu=-4, tau=(1 / 2**2))]
# pos: the first pos is just the observation
pos_states = [pm.Normal("p0", mu=observations[0], tau=tau_system_noise)]

for i in range(1, len(observations)):
    new_vel = pm.Normal("v" + str(i),
                        mu=vel_states[-1] + damping * dt,
                        tau=tau_system_noise)
    vel_states.append(new_vel)
    pos_states.append(
        pm.Normal("s" + str(i),
                  mu=pos_states[-1] + new_vel * dt,
                  tau=tau_system_noise)
    )

# we assume some observation noise
tau_observation_noise = (1 / 0.5**2)
obs = pm.Normal("obs", mu=pos_states, tau=tau_observation_noise, value=observations, observed=True)

これは私がサンプリングを実行する方法です：

mcmc = pm.MCMC([damping, obs, vel_states, pos_states])
mcmc.sample(50000, 25000)
pm.Matplot.plot(mcmc.get_node("damping"))
damping_samples = mcmc.trace("damping")[:]
print "damping -- mean:%f; std:%f" % (mean(damping_samples), std(damping_samples))
print "real damping -- %f" % true_damping

の値dampingは事前確率によって支配されます。ユニフォームなどの前を変えてもそのままです。

私は何を間違っていますか？別のレイヤーがあるだけで、前の例とほとんど同じです。

この問題の完全な IPython ノートブックは、http: //nbviewer.ipython.org/github/sotte/random_stuff/blob/master/PyMC%20-%20HMM%20Dynamic%20System.ipynbで入手できます。

[編集: いくつかの説明とサンプリングのためのコード。]

[EDIT2: @Chris の回答は役に立ちませんでした。AdaptiveMetropolis*_states はモデルの一部ではないため、使用できませんでした。]

4

4 に答える 4

1

モデルの構造を知っていると仮定すると (システム同定ではなくパラメーター推定を行っている)、未知の減衰、初期位置、初期速度をパラメーターとして、位置の配列、観測値を使用して、PyMC モデルを回帰として構築できます。 .

つまり、点-質量系を表すクラス PM を使用すると、次のようになります。

pm = PM(true_damping)

positions, velocities = pm.integrate(true_pos, true_vel, N, dt)

# Assume little system noise
std_system_noise = 0.05
tau_system_noise = 1.0/std_system_noise**2

# Treat the real positions as observations
observations = positions + np.random.randn(N,)*std_system_noise

# Damping is modelled with a Uniform prior
damping = mc.Uniform("damping", lower=-4.0, upper=4.0, value=-0.5)

# Initial position & velocity unknown -> assume Uniform priors
init_pos = mc.Uniform("init_pos", lower=-1.0, upper=1.0, value=0.5)
init_vel = mc.Uniform("init_vel", lower=0.0, upper=2.0, value=1.5)

@mc.deterministic
def det_pos(d=damping, pi=init_pos, vi=init_vel):
    # Apply damping, init_pos and init_vel estimates and integrate
    pm.damping = d.item()
    pos, vel = pm.integrate(pi, vi, N, dt)
    return pos

# Standard deviation is modelled with a Uniform prior
std_pos = mc.Uniform("std", lower=0.0, upper=1.0, value=0.5)

@mc.deterministic
def det_prec_pos(s=std_pos):
    # Precision, based on standard deviation
    return 1.0/s**2

# The observations are based on the estimated positions and precision
obs_pos = mc.Normal("obs", mu=det_pos, tau=det_prec_pos, value=observations, observed=True)

# Create the model and sample
model = mc.Model([damping, init_pos, init_vel, det_prec_pos, obs_pos])
mcmc = mc.MCMC(model)
mcmc.sample(50000, 25000)

完全なリストはこちら: https://gist.github.com/stukeyr/7762371

N を増やして dt を減らすと、推定値が大幅に改善されます。

于 2013-12-03T01:46:03.413 に答える
1

もう一度見てみると、モデルにはいくつかの問題があります。何よりもまず、すべての PyMC オブジェクトをモデルに追加していません。のみを追加しまし[damping, obs]た。すべての PyMC ノードをモデルに渡す必要があります。

Modelまた、 と の両方を呼び出す必要はないことに注意してくださいMCMC。これで問題ありません:

model = pm.MCMC([damping, obs, vel_states, pos_states])

PyMC の最適なワークフローは、モデルを実行中のロジックとは別のファイルに保持することです。そうすれば、モデルをインポートしてに渡すことができますMCMC：

import my_model

model = pm.MCMC(my_model)

または、モデルを関数として記述し、locals(またはvars) を返し、その関数を の引数として呼び出すこともできますMCMC。例えば：

def generate_model():

    # put your model definition here

    return locals()

model = pm.MCMC(generate_model())
于 2013-12-02T15:15:29.260 に答える
0

不合理とはどういう意味ですか？それらは事前に縮小されていますか？ダンピングの分散はかなり狭いようです。より拡散した事前設定を与えるとどうなるでしょうか?

また、AdaptiveMetropolis*_states 配列でサンプラーを使用してみてください。

my_model.use_step_method(AdaptiveMetropolis, my_model.vel_states)

相関変数のほうがよく混ざる場合があります。

于 2013-11-30T20:52:09.230 に答える