以log(X)为状态变量的Python扩展卡尔曼滤波实现疑问

问题描述 #

我正在尝试实现一个扩展卡尔曼滤波（Extended Kalman Filter），希望将其中一个状态变量（如x₃）建模为log(x₃)以确保其非负性。现有普通扩展卡尔曼滤波代码如下：

# Loop over time
for t in range(T):
    
    # Predict state vector
    Xt_1[:, [t]] = muP + thetaP @ Xt[:, [t]]
    
    # MSE matrix state equation
    Pt_1[:, :, t] = thetaP @ Pt[:, :, t] @ thetaP.T + Sigmae
    
    # Measurement equation
    y_hat = A + B.T @ Xt_1[:, [t]]
    # fitting errors
    errors = yields[[t], :] - y_hat
    
    # estimatge Jacobian
    H = compute_jacobian(Xt_1, t)
    # compute F
    F = np.matrix(H @ Pt_1[:, :, t] @ H.T) + Sigmaz

    # Kalman gain
    K = Pt_1[:, :, t] @ B @ np.linalg.inv(F)
    # state vector
    Xt[:, [t + 1]] = Xt_1[:, [t]] + K @ errors
    # MSE matrix state equation
    Pt[:, :, t + 1] = (np.eye(3) - K @ B.T) @ Pt_1[:, :, t]

若将状态向量Xt设为(x₁, x₂, log(x₃))而非原有的(x₁, x₂, x₃)，除修改雅可比矩阵和设置log(x₃)初始值外，是否需要修改预测和更新阶段的状态方程？

解答 #

是的，除了修改雅可比矩阵和初始值，你还需要对预测阶段的状态方程、测量方程，以及更新阶段的相关计算做调整，具体如下：

1. 预测阶段的状态方程调整 #

原代码的预测基于线性状态转移Xt_1 = muP + thetaP @ Xt，但现在状态变量第三项是log(x₃)，需根据原x₃的转移逻辑调整：

• 如果原状态转移中x₃是线性更新（比如x₃(t+1) = a*x₃(t) + b + 噪声），转换到log(x₃)后，状态转移会变成非线性。比如原转移是x₃' = c1*x₁ + c2*x₂ + c3*x₃ + mu3 + e，替换为z3 = log(x₃)后，方程变为exp(z3') = c1*x₁ + c2*x₂ + c3*exp(z3) + mu3 + e，此时不能直接用原来的thetaP矩阵，需要重新推导非线性状态转移函数，并用状态转移的雅可比矩阵代替thetaP计算预测协方差Pt_1。
• 如果原状态转移本身就是针对log(x₃)设计的线性模型，可保留类似形式，但要确保thetaP的第三行/列对应log(x₃)的转移系数。

2. 测量方程的调整 #

原测量方程直接用包含log(x₃)的Xt_1计算y_hat是错误的，因为实际测量模型中需要的是x₃而非其对数，必须先转换回原变量：

# 从预测状态中提取变量，转换log(x₃)为x₃
x1_pred = Xt_1[0, t]
x2_pred = Xt_1[1, t]
z3_pred = Xt_1[2, t]
x3_pred = np.exp(z3_pred)
# 重构用于测量的状态向量
X_meas_pred = np.array([[x1_pred], [x2_pred], [x3_pred]])
# 计算预测测量值
y_hat = A + B.T @ X_meas_pred

3. 雅可比矩阵计算逻辑修改 #

原compute_jacobian是针对(x₁,x₂,x₃)计算测量方程的雅可比矩阵，现在要适配新状态变量(x₁,x₂,log(x₃))：
假设测量方程对x₃的偏导是dy/dx3，那么对log(x₃)的偏导就是dy/dx3 * exp(z3)（因为dx3/d(log(x3)) = x3 = exp(z3)），所以compute_jacobian函数需要重新推导每个元素的偏导数。

4. 更新阶段的协方差调整（可选但重要） #

原代码的协方差更新Pt[:, :, t+1] = (np.eye(3) - K @ B.T) @ Pt_1[:, :, t]是线性场景下的简化形式。如果状态转移变为非线性，标准EKF的协方差更新需要基于状态转移的雅可比矩阵（比如用Phi表示）来计算预测协方差，更新阶段的公式也需要对应调整；如果状态转移对新变量仍为线性，这部分可保留。

总结 #

核心修改点：

• 预测阶段：根据原x₃的转移逻辑，重新推导针对log(x₃)的状态转移函数（可能从线性变非线性）
• 测量阶段：将log(x₃)转换回x₃后再代入测量方程计算y_hat
• 雅可比矩阵：重新计算针对新状态变量的测量雅可比和状态转移雅可比

内容的提问来源于stack exchange，提问作者Jessica F.