隐马尔可夫模型的前向-后向算法中的数值不稳定性

Question

我正在实现隐马尔可夫模型的前向算法（算法见下文）。为了防止上溢/下溢，我使用对数概率，并使用对数和表达式技巧来计算每个前向系数。

我绘制了计算出的前向系数，并将其与我用来模拟数据的状态进行比较。如下图所示，总体形状看起来是正确的，因为前向系数在与状态相同的位置出现峰值。问题是前向系数是概率，因此它们的对数永远不应该超过 0，但是在下图中，我看到存在逐渐漂移，并且对数系数明显超过零，我怀疑这是由于累积的数值误差造成的。 （注意，在我的符号中，g_j(z_j) 表示状态 z_j=1 或 2 时的前向系数在 j 时刻的对数）。

我已经使用了 log-sum-exp 技巧，所以我想知道我还能做些什么来解决这个问题？ （防止对数概率超过 0，并消除这种逐渐向上的漂移）。

我的代码的相关部分如下：

    def log_sum_exp(self, sequence):
        '''
        Returns np.log( np.sum(sequence) ) without under/overflow.
        '''
        sequence = np.array(sequence)
        if np.abs(np.min(sequence)) > np.abs(np.max(sequence)):
            b = np.min(sequence)
        else:
            b = np.max(sequence)

        return b + np.log(np.sum(np.exp(sequence-b)))

    def g_j_z(self, j, z_j):
        '''
        Returns g_j(z_j).
        j: (int) time index. zero-indexed 0, 1, 2, ... n-1
        z_j: (int) state index. zero-indexed. 0, 1, 2, ... K-1
        '''
        
        if j == 0:
            return np.log(self.p_init[z_j]) + self.log_distributions[z_j](self.pre_x+[self.x[0]], self.pre_exog+[self.exog[0]])

        if (j, z_j) in self.g_cache:
            return self.g_cache[(j, z_j)]

        temp = []
        for state in range(self.K):
            temp.append(
                self.g_j_z(j-1, state) + np.log(self.p_transition[state][z_j]) 
            )

        self.g_cache[(j, z_j)] = self.log_sum_exp(temp) + self.log_distributions[z_j](self.pre_x+self.x[0:j+1], self.pre_exog+self.exog[0:j+1])  

        return self.g_cache[(j, z_j)]

变量解释：

self.g_cache

是一个字典，将元组

(j, z_j)

（时间和状态）映射到对数系数 g_j(z_j)。这是为了避免重复计算。

self.p_init

是一个列表。

self.p_init[i]

包含处于状态

i

的初始概率。

self.p_transition

是一个矩阵。

self.p_transition[i][j]

包含从状态

i

转换到状态

j

的概率。

self.log_distributions

是函数列表。

self.log_distributions[i]

是状态

i

的对数概率分布，它是一个以观测历史和外生变量作为输入的函数，并返回最新观测的对数概率。例如，对于一个AR-1进程，日志分发的实现如下

def log_pdf1(x, exog, params=regime1_params):
    '''
    x: list of all history of x up to current point
    exog: list of all history of exogenous variable up to current point
    '''
    # AR1 process with exogenous mean
    alpha, dt, sigma = params[0], params[1], params[2]
    mu = x[-2] + alpha*(exog[-1] - x[-2])*dt
    std = sigma*np.sqrt(dt)
    return norm.logpdf(x[-1], loc=mu, scale=std)

我正在实现的算法在这里给出：

但是，我使用 log-sum-exp 技巧来计算系数的对数，以避免上溢/下溢：

非常感谢您的帮助！

Answer 1

前向系数对数的正漂移可能是由于随着时间的推移累积数值误差造成的，这是重复应用 log-sum-exp 技巧时的常见问题。由于对数概率应保持非正值，因此任何高于零的增加都表明每次迭代都会产生小的不准确性。为了缓解这种情况，请考虑在每个时间步标准化前向系数。例如，您可以通过减去每一步的最大对数系数来重新居中对数概率。这种标准化在不改变序列的情况下保持相对值，有助于控制数值误差并将值保持在预期的 范围内。