拟合切换线性动态系统（SLDS）时遭遇运行时警告问题

原始代码 #

import matplotlib.pyplot as plt
# 注意：原代码缺失numpy和ssm导入，运行会报错
import numpy as np
import ssm

# 时间序列长度
time = 100

# 两组动力学参数
A1 = np.array([[0.9, 0], [0, 0.9]])  # A1对应稳定节点
A2 = np.array([[0.9, 0.5], [-0.5, 0.9]])  # A2对应稳定螺旋

# 初始化状态和观测值
x = np.zeros((time, 2)) 
y = np.zeros((time, 2))

# 生成数据
for t in range(1, time):
    if t < 50:
        x[t] = A1 @ x[t-1] + np.random.normal(0, 0.1, 2) 
    else:
        x[t] = A2 @ x[t-1] + np.random.normal(0, 0.1, 2)
    y[t] = x[t] + np.random.normal(0, 0.1, 2)  # 带观测噪声的观测值

# 绘制生成的数据
plt.plot(y[:,0], label='维度1') # y[:,0]表示选取y的所有行第一列数据，Python中:代表从起始到结束的完整切片
plt.plot(y[:,1], label='维度2')
plt.legend()
plt.show()

# 设置SLDS模型参数
time_bins = 100    # 时间步数
n_disc_states = 2       # 离散状态数
latent_dim = 2       # 隐变量维度
emissions_dim = 2      # 观测维度

# 初始化SLDS模型
model = ssm.SLDS(emissions_dim, n_disc_states, latent_dim, emissions="gaussian") 

# 使用BBVI结合均值场后验方法拟合模型
# 注意：原代码使用未定义的y_scaled，需替换为y或先做归一化
elbos, posterior = model.fit(y_scaled, method="bbvi", variational_posterior="mf", num_iters=1000, init_method="kmeans")

遇到的警告 #

RuntimeWarning: divide by zero encountered in log
return f_raw(*args, **kwargs)
RuntimeWarning: divide by zero encountered in divide
lls = -0.5 * np.log(2 * np.pi * etas) - 0.5 * (data[:, None, :] - mus)**2 / etas
RuntimeWarning: invalid value encountered in subtract
lls = -0.5 * np.log(2 * np.pi * etas) - 0.5 * (data[:, None, :] - mus)**2 / etas

问题根源 #

1. 噪声参数初始化异常：SSM库的高斯发射模型中，噪声参数etas被初始化为极小值甚至0，导致np.log(etas)触发除零警告，后续除法、减法运算产生无效值。
2. 代码缺失依赖与未定义变量：原代码未导入numpy和ssm库，且使用了未定义的y_scaled变量，这些都会导致运行异常。

修复方案 #

• 补全依赖导入：在代码开头添加import numpy as np和import ssm，确保依赖齐全。
• 修复未定义变量：
  • 若无需归一化，直接将model.fit(y_scaled, ...)改为model.fit(y, ...)；
  • 若需归一化数据，先执行归一化处理：

    from sklearn.preprocessing import StandardScaler
    scaler = StandardScaler()
    y_scaled = scaler.fit_transform(y)
    

• 限制噪声参数范围：初始化模型后，手动约束etas的最小值，避免趋近于0：

  # 初始化SLDS模型后添加
  model.emissions.etas = np.clip(model.emissions.etas, 1e-6, None)
  

• 更换初始化方法：kmeans初始化可能导致参数分布极端，换成meanfield或random提高稳定性：

  elbos, posterior = model.fit(y, method="bbvi", variational_posterior="mf", num_iters=1000, init_method="meanfield")
  

• 调整BBVI学习率：添加learning_rate参数，避免参数更新幅度过大导致噪声参数异常：

  elbos, posterior = model.fit(y, method="bbvi", variational_posterior="mf", num_iters=1000, init_method="kmeans", learning_rate=1e-3)
  

内容的提问来源于stack exchange，提问作者sssj