PyMC: иерархические Hidden Markov Model

Question

Это является продолжением на PyMC: Parameter estimation in a Markov system

У меня есть система, которая определяется его положением и скоростью на каждом временном шаге. Поведение системы определяется как:

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

[EDIT:. Некоторые пояснения & кода для отбора проб]

[EDIT2: @ Крис ответ не помог. Я не мог использовать AdaptiveMetropolis, так как * _states не являются частью модели.]

Answer 1

Что вы подразумеваете под необоснованным? Они урезаны по отношению к предыдущему? У демпфирования, похоже, довольно сложная дисперсия - что, если вы дадите ему более диффузное ранее?

Кроме того, вы можете попробовать использовать AdaptiveMetropolis пробник на * _states массивы:

my_model.use_step_method(AdaptiveMetropolis, my_model.vel_states) 

Это иногда смешивается лучше для коррелированных переменных, так как они, вероятно, являются.

Answer 2

Есть пара проблем с моделью, снова глядя на нее. Прежде всего, вы не добавили все ваши объекты 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()) 

Answer 3

Предполагая, что вы знаете структуру вашей модели - вы делаете оценку параметров, а не идентификации системы - вы можете построить модель 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/stuckeyr/7762371

Увеличение N и уменьшение DT улучшит ваши оценки заметно.

Answer 4

Я думаю, что ваш первоначальный подход прекрасен и должен работать, за исключением того, что переменная «obs» не включена в список узлов, поставляемых в MCMC (см. Раздел [10] в вашем ноутбуке). После включения этой переменной решатель MCMC работает отлично и выполняет принудительные зависимости, определенные вашей моделью. Я хотел бы повторить точку зрения Криса, что лучше всего определить модель в другом файле или под функцией, чтобы избежать таких ошибок.

Причина, по которой вы не получаете правильных результатов, заключается в том, что ваши приоритеты были выбраны произвольно, а в некоторых случаях значения таковы, что модели очень сложно правильно смешивать, чтобы сходиться. Ваша игрушечная проблема пытается оценить величину демпфирования, так что позиции сходятся к вектору наблюдаемых положений. Для этого ваша модель должна иметь гибкость в выборе значений скорости и демпфирования в широком диапазоне, чтобы стохастические ошибки в позиции/скорости можно было исправить при переходе от одного временного шага к следующему. В противном случае, в результате вашей схемы интеграции Эйлера, ошибки просто продолжают распространяться. Я думаю, что Крис упомянул то же самое, когда он предложил выбрать более диффузный ранее.

Предлагаю поиграть с значениями tau для каждой из обычных переменных. Например, я изменил следующие значения:

damping = pm.Normal("damping", mu=0, tau=1/20.**2) # was tau=1/2.**2 

new_vel = pm.Normal("v" + str(i), 
        mu=vel_states[-1] + damping * dt, 
        tau=(1/2.**2)) # was tau=tau_system_noise=(1/0.5**2) 

tau_observation_noise = (1/0.005**2) # was 1/0.5**2 

Вы можете увидеть измененный файл here.

Участки внизу показывают, что позиции действительно сходятся. Скорости повсюду. Оценочное среднее значение демпфирования составляет 6,9, что сильно отличается от -1,5. Возможно, вы сможете добиться более высоких оценок, выбирая подходящие значения для priors.