Численная нестабильность в алгоритме вперед-назад для скрытых марковских моделей ⇐ Python

[Anonymous]

Я реализую прямой алгоритм для скрытых марковских моделей (алгоритм см. ниже). Чтобы предотвратить переполнение/недополнение, я вместо этого работаю с логарифмическими вероятностями и использую прием log-sum-exp для вычисления каждого прямого коэффициента.
Я построил вычисленный прямой коэффициент и сравнил его с состояния, которые я использовал для моделирования своих данных. Как показано на рисунке ниже, общая форма выглядит правильной, поскольку прямой коэффициент резко возрастает в тех же местах, что и состояния. Проблема в том, что прямые коэффициенты являются вероятностями, поэтому их журналы никогда не должны превышать 0, однако на изображениях ниже я вижу постепенный дрейф, и логарифмические коэффициенты явно превышают ноль, что, как я подозреваю, связано с накопленными числовыми ошибками.(Обратите внимание, в моих обозначениях g_j(z_j) обозначает логарифм прямого коэффициента в момент времени j для состояния z_j=1 или 2).


Я уже использовал трюк 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.

Код: Выделить всё

self.p_transition

— это матрица. self.p_transition[j] содержит вероятность перехода из состояния i в состояние j.

Код: Выделить всё

self.log_distributions

— список функций. self.log_distributions — это лог-распределение вероятностей для состояния 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)

Алгоритм, который я реализую, приведен здесь:
[img]https://i. sstatic.net/8MHvufYT.png[/img]

Однако вместо этого я вычисляю журнал коэффициентов, используя трюк log-sum-exp, чтобы избежать переполнения/недополнения:


Большое спасибо за помощь!

Подробнее здесь: https://stackoverflow.com/questions/791 ... kov-models