当前位置：首页 > news >正文

输入驱动HMM实战：用SSM构建GLM-HMM模型，揭示外部 covariates如何影响状态转移

news 2026/8/10 20:40:16

输入驱动HMM实战：用SSM构建GLM-HMM模型，揭示外部 covariates如何影响状态转移

【免费下载链接】ssmBayesian learning and inference for state space models项目地址: https://gitcode.com/gh_mirrors/ssm38/ssm

SSM（State Space Models）是一个强大的贝叶斯学习与推理工具包，专注于状态空间模型的构建与分析。本文将带你通过实战案例，学习如何使用SSM构建GLM-HMM（广义线性模型-隐马尔可夫模型），深入理解外部 covariates如何影响系统的状态转移过程。

什么是GLM-HMM？

GLM-HMM是一种融合了广义线性模型（GLM）和隐马尔可夫模型（HMM）的混合模型。与传统HMM不同，GLM-HMM假设观测概率依赖于外部输入 covariates，而状态转移概率保持平稳。这种结构特别适合分析外部因素如何影响系统的动态行为，例如动物在决策任务中的行为模式。

在GLM-HMM中，观测概率的计算方式如下：

对于类别c ≠ C： $$ \Pr(y_t = c \mid z_{t} = k, u_t, w_{kc}) = \frac{\exp{w_{kc}^\mathsf{T} u_t}}{1+\sum_{c'=1}^{C-1} \exp{w_{kc'}^\mathsf{T} u_t}} $$

对于类别c = C： $$ \Pr(y_t = c \mid z_{t} = k, u_t, w_{kc}) = \frac{1}{1+\sum_{c'=1}^{C-1} \exp{w_{kc'}^\mathsf{T} u_t}} $$

其中，$u_{t}$ 是输入协变量，$w_{kc}$ 是与状态k和类别c相关的权重参数。

快速入门：用SSM构建GLM-HMM模型

环境准备

首先，确保你已经安装了SSM库。如果还没有安装，可以通过以下命令获取项目代码并安装：

git clone https://gitcode.com/gh_mirrors/ssm38/ssm cd ssm pip install -e .

初始化GLM-HMM模型

使用SSM构建GLM-HMM模型非常简单。以下代码展示了如何创建一个具有输入驱动观测的HMM：

import ssm # 设置模型参数 num_states = 3 # 离散状态数量 obs_dim = 1 # 观测维度 num_categories = 2 # 输出类别数量 input_dim = 2 # 输入维度 # 创建GLM-HMM模型 glmhmm = ssm.HMM(num_states, obs_dim, input_dim, observations="input_driven_obs", observation_kwargs=dict(C=num_categories), transitions="standard")

实战案例：模拟动物决策行为

生成模型参数

我们模拟一个动物在2AFC（二选一）决策任务中的行为。设置3个状态，每个状态具有不同的GLM权重和转移概率：

# 设置生成模型的权重和转移矩阵 gen_weights = np.array([[[6, 1]], [[2, -3]], [[2, 3]]]) gen_log_trans_mat = np.log(np.array([[[0.98, 0.01, 0.01], [0.05, 0.92, 0.03], [0.03, 0.03, 0.94]]])) # 应用参数到模型 true_glmhmm.observations.params = gen_weights true_glmhmm.transitions.params = gen_log_trans_mat

模拟输入序列和观测数据

我们生成模拟的外部输入序列，包括刺激值和偏置项，然后使用模型生成观测数据：

# 创建外部输入序列 num_sess = 20 # 会话数量 num_trials_per_sess = 100 # 每个会话的试次数 inpts = np.ones((num_sess, num_trials_per_sess, input_dim)) stim_vals = [-1, -0.5, -0.25, -0.125, -0.0625, 0, 0.0625, 0.125, 0.25, 0.5, 1] inpts[:,:,0] = np.random.choice(stim_vals, (num_sess, num_trials_per_sess)) inpts = list(inpts) # 生成状态和观测数据 true_latents, true_choices = [], [] for sess in range(num_sess): true_z, true_y = true_glmhmm.sample(num_trials_per_sess, input=inpts[sess]) true_latents.append(true_z) true_choices.append(true_y)

模型拟合与参数恢复

最大似然估计（MLE）

使用EM算法拟合模型：

# 创建新的GLM-HMM模型 new_glmhmm = ssm.HMM(num_states, obs_dim, input_dim, observations="input_driven_obs", observation_kwargs=dict(C=num_categories), transitions="standard") # 拟合模型 N_iters = 200 fit_ll = new_glmhmm.fit(true_choices, inputs=inpts, method="em", num_iters=N_iters, tolerance=10**-4)

状态 permutation 与参数比较

由于HMM状态具有可交换性，我们需要对状态进行permutation以与真实模型比较：

from ssm.util import find_permutation # 对状态进行permutation new_glmhmm.permute(find_permutation(true_latents[0], new_glmhmm.most_likely_states(true_choices[0], input=inpts[0]))) # 比较恢复的权重和转移矩阵与真实参数 recovered_weights = new_glmhmm.observations.params recovered_trans_mat = np.exp(new_glmhmm.transitions.log_Ps)

高级技巧：最大后验估计（MAP）与交叉验证

使用MAP估计处理小数据问题

在数据量有限时，MAP估计可以通过引入先验知识提高模型泛化能力：

# 创建带先验的GLM-HMM模型 prior_sigma = 2 prior_alpha = 2 map_glmhmm = ssm.HMM(num_states, obs_dim, input_dim, observations="input_driven_obs", observation_kwargs=dict(C=num_categories, prior_sigma=prior_sigma), transitions="sticky", transition_kwargs=dict(alpha=prior_alpha, kappa=0)) # 用MAP估计拟合模型 _ = map_glmhmm.fit(true_choices, inputs=inpts, method="em", num_iters=N_iters, tolerance=10**-4)

交叉验证选择最佳模型

通过交叉验证比较MLE和MAP估计的性能：

# 生成测试数据 num_test_sess = 10 test_inpts = np.ones((num_test_sess, num_trials_per_sess, input_dim)) test_inpts[:,:,0] = np.random.choice(stim_vals, (num_test_sess, num_trials_per_sess)) test_inpts = list(test_inpts) # 比较模型在测试数据上的似然度 mle_test_ll = new_glmhmm.log_likelihood(test_choices, inputs=test_inpts) map_test_ll = map_glmhmm.log_likelihood(test_choices, inputs=test_inpts)

扩展应用：多项GLM-HMM

SSM同样支持多类别输出的GLM-HMM模型。以下是一个具有3个输出类别的示例：

# 设置多项GLM-HMM参数 num_states = 4 num_categories = 3 input_dim = 2 # 创建多项GLM-HMM模型 multinom_glmhmm = ssm.HMM(num_states, obs_dim, input_dim, observations="input_driven_obs", observation_kwargs=dict(C=num_categories), transitions="standard")