在PyMC3中再现哈密顿1989马尔可夫切换模型

import pymc3 as pm
import numpy as np
import matplotlib.pyplot as plt

# 定义状态转移矩阵
transition_matrix = np.array([[0.9, 0.1], [0.2, 0.8]])

# 定义初始状态分布
initial_state = np.array([0.5, 0.5])

# 定义观测分布
observation_distribution = np.array([[0.6, 0.4], [0.3, 0.7]])

# 定义模型
with pm.Model() as model:
    # 定义状态序列
    states = pm.Categorical('states', p=initial_state, shape=len(data))
    
    # 定义观测序列
    observations = pm.Categorical('observations', p=observation_distribution[states], observed=data)
    
    # 定义状态转移
    for i in range(len(data) - 1):
        pm.MarkovSwitching('state_transition_{}'.format(i), p=transition_matrix[states[i]], observed=states[i+1])

with model:
    trace = pm.sample(1000, tune=1000)

pm.plot_posterior(trace, var_names=['transition_matrix', 'initial_state', 'observation_distribution'])
plt.show()

mode_states = np.argmax(np.bincount(trace['states'], minlength=len(data)))
print('Mode states:', mode_states)