스칼라 칼만 자료동화

1. 시간에 따라 바뀌는 모형에 관측 결합하기

자료동화는 새 관측이 들어올 때 상태 추정값을 갱신합니다. 스칼라 상태의 동역학을 $x_{n+1}=a x_n+w_n$, 관측을 $y_{n+1}=x_{n+1}+v_{n+1}$로 둡니다. 과정 잡음과 관측 잡음은 서로 독립인 평균 0의 가우스 잡음이며 분산은 $q$와 $r$입니다. 초기 상태도 가우스 분포이면 이 모형의 조건부 평균과 분산을 칼만 필터로 정확히 계산할 수 있습니다. Kalman의 선형 필터링 및 예측 연구

2. 다음 상태와 불확실성 예측하기

이전 분석의 평균이 $m_n$, 분산이 $P_n$이면 다음 예측은

입니다. 상태에 $a$를 곱하면 분산은 $a^2$배가 됩니다. 독립적인 과정 잡음의 분산 $q$를 더합니다. 다음 관측이 없어도 예측 단계는 실행할 수 있습니다.

3. 관측과 예측의 차이로 갱신하기

새 관측과 예측의 차이 $y-m^-$를 혁신이라 하며, 모형이 예상하는 분산은 $P^-+r$입니다. 칼만 이득 $K=P^-/(P^-+r)$로 분석을 갱신합니다.

관측 잡음이 크면 이득은 작아집니다. 관측이 빠졌다면 $m^-,P^-$를 그대로 다음 단계에 전달합니다. 과정 잡음으로 관측 공백 중 분산이 커질 수 있습니다. 모형의 수축이 강하고 과정 잡음이 작으면 예측 분산이 줄 수도 있습니다.

4. 첫 예측과 분석 손으로 계산하기

$m_0=0$, $P_0=1$, $a=0.9$, $q=0.2$, $r=0.25$로 둡니다. 첫 예측 분산은 $P^-=1.01$입니다. $y_1=1$이면 혁신은 1이고 $K=1.01/1.26\approx0.8016$입니다. 분석 평균은 약 0.8016, 분산은 약 0.2004입니다. 다음 관측이 없으면 이 값에서 예측만 계산합니다.

5. 관측 하나를 뺀 필터를 Python으로 실행하기

from math import isclose

a, q, r = .9, .2, .25
m, p = 0., 1.
observations = [1., None, .2]
for t, y in enumerate(observations, 1):
    m_pred, p_pred = a*m, a*a*p+q
    if y is None:
        m, p = m_pred, p_pred
        gain = 0.
    else:
        gain = p_pred/(p_pred+r)
        m = m_pred+gain*(y-m_pred)
        p = p_pred*r/(p_pred+r)
        assert p < p_pred
    if t == 1:
        assert isclose(m, 1.01/1.26)
    print(t, y, f"gain={gain:.4f}, mean={m:.4f}, variance={p:.4f}")

6. 관측과 모형이 어긋나는 부분 찾기

관측이 들어올 때마다 혁신과 예측 분산을 기록합니다. 혁신의 부호가 한쪽으로 계속 치우치면 편향이나 상태 모형의 부족을 의심할 수 있습니다. 조건별로 정규화한 혁신을 비교하고 체계적인 패턴이 있는지 확인합니다. 필터는 현재까지의 관측으로 각 시각의 상태를 추정합니다. 평활화는 이후 관측도 사용해 앞선 상태 추정을 수정할 수 있습니다.