[시계열의 상태공간 모델] Kalman Filter

State Space Models for Time Series, 시계열의 상태공간 모델

 

일반적으로 사용되는 State space models는 다음과 같다.

• The Kalman filter applied to a linear Gaussian model (선형 가우스 모델에 적용된 칼만 필터)
• Hidden Markov Models (은닉 마르코프 모형)
• Bayesian structural time series (베이즈 구조적 시계열)

각 방법에 대한 모델은 이미 잘 구현되어 있으므로, 각 모델에 대한 수학적인 직관력을 키우고 어떤 데이터가 적합한지에 대한 내용을 살펴볼 것이다. 실제 코드 예제 또한 다룰 예정이다.

 

여기서 '관측한 것'과 '관측이 만들어낸 상태'를 구분할 것이다. 이때, 관측에 기반한 underlying state을 추정하는 작업은 다음 단계로 나타낼 수 있다.

• Filtering(필터링): 시간 t의 상태에 대한 추정을 update하는 것에 시간 t에서의 measurement를 사용
  • 상태에 대한 추정을 update할 때 가장 최근 정보에 가중치를 주는 방법
• Forecasting(예측): 시간 t의 예상되는 상태에 대한 예측을 하는 데에 시간 t-1에서의 measurement를 사용
  • 미래에 대한 정보 없이 미래를 예측
• Smoothing(평활화): 시간 t의 '참 상태' 추정에 시간 t 및 t의 전후를 사용
  • 특정 시간에서의 best 상태 추정에 미래와 과거 정보 모두 사용

이 작업들은 유사해보이지만, 차이를 구분하는 것이 중요하다.

 

---

State Space Models(상태공간) 모델의 장단점

장점

• deterministic하거나 stochastic한 적용을 할 수 있다.
• continuous하거나 discrete하게 sampling된 데이터에도 적용될 수 있다.
• 노이즈가 낀 데이터 자체보다, 노이즈한 데이터의 생성 과정과 상태를 모델링할 수 있게 해준다.
• process의 생성을 설명할 수 있는 인과관계에 대한 모델을 모델링 단계에서 주입한다.
• 이미 익숙한 시스템의 dynamics를 디테일하게 파악할 수 있다.

상태공간 모델은 시간에 따른 상관계수, 파라미터의 변화를 허용한다. 즉 시간에 따라 변화하는 행동을 허용한다는 것이다. 또한 데이터에 정상성이라는 조건을 부과하지 않는다.

 

 

단점

• 모델이 너무 flexible해서 수많은 파라미터가 설정되고, 수많은 종류의 모델이 존재할 수 있다. (즉 특정 상태공간 모델을 구성하는 모든 properties가 심도 있게 연구되지 않음, 시계열의 상태공간 모델에 대한 선행연구가 부족)
• 수많은 파라미터로 구성된 상태공간 모델은 계산 측면에서 부담스럽다. (또한 일부 모델은 과적합에 취약. 데이터가 많지 않을 수록 더 그럴 수 있다)

 

---

The Kalman filter applied to a linear Gaussian model (칼만 필터)

칼만 필터는 시계열로부터 새로운 정보를 결합하는 데에 사용되고, underlying state을 추정하기 위해 이전의 정보들에 새로운 정보를 결합하는 좋은 방법이다. (칼만 필터가 최초로 사용된 사례로 아폴로 11호 미션이 있다. 기내 계산 자원만으로는 위치추정이 불가하기에, 계산이 비교적 간단하고, 미래 예측이나 현재를 추정할 때에 과거 데이터의 저장이 필요 없기 때문)

 

 

그럼 공식을 한 번 이해해보자.

먼저 state와 observation이 다음의 dynamics를 가진다고 상정하는 Linear Gaussian model로 시작한다.

즉 시간 t에서의 상태(x_t)는 이전 시간의 상태 항인 F * x_(t-1), 외부 압력에 대한 항인 B * u_t, 확률적 항인 w_t에 대한 함수로 표현할 수 있다.

이와 비슷하게 시간 t에서의 측정(y_t)은 시간 t의 상태, 확률적 오차항, 측정 오차항의 함수다.

 

만일 x_t가 우주선의 실제 위치를 나타내고, y_t가 센서로 측정한 위치라고 생각해보자.

이때 v_t는 센서장치의 측정 오차를 의미한다. 그러면 칼만 필터에 적용 가능한 기본 공식은 다음과 같다.

시간 t에서 새로운 정보가 주어졌을 때 추정을 update하는 방법을 보여주는데, 이는 필터링 단계를 말한다. 즉 시간 t의 상태에 대한 추정을 update하는 데에 시간 t의 측정을 사용하는 것이다.

(여기서 알아두어야 할 것은, y_t만 관측이 가능하며 정확한 상태는 모른 채 상태를 추론했다는 상황을 가정한다)

그리고 위 식의 K_t는 새로운 정보(y_t)와 과거의 정보( (x hat)_(t-1) ) 사이의 추정에 대한 균형을 맞추기 위해 사용된 것임을 알 수 있다.

 

모델의 동작 방식을 좀 더 깊게 살펴보자.

• state의 공분산 추정은 P_t로 표현 (state가 단변량인지, 다변량인지에 따라 스칼라 혹은 행렬이 될 수 있는데, 대체로 다변량이 일반적)
• (P-)_t는 시간 t에서의 측정이 고려되기 전의 추정
• 측정 오차의 분산은 R (즉 v_t에 대한 분산. 일반적으로 꽤 잘 정의되어 있고, 마찬가지로 측정의 차원에 따라 스칼라 혹은 공분산 행렬이 될 수 있음)
• w_t의 적절한 값인 Q는 R과는 달리, 명확하지 않으며 모델링 시 조정이 필요

그러면 시간 0에서 이미 알고 있는 x와 추정된 P로 구성된 process에서 시작한다. 그 뒤, 시간 0에서 시간 뒤로 움직이며 prediction(예측)과 update(갱신)을 반복한다.

먼저 prediction 단계가 발생한 뒤, update / filtering 단계가 실행되고, 이 순서대로 계속 반복한다.

 

Prediction

 

Filtering

 

 

여기서 K_t는 the Kalman gain(칼만 이득)으로, 다음을 뜻한다.

 

 

설명하자면

(1) 측정 y_t에 대한 내용을 알지못한 채, 시간 t에서의 예측에 수행된 계산 단계(예측)와

(2) y_t를 알게된 이후 시간 t에서 수행된 단계

가 있다는 것이다.

 

이 과정을 진행하기 위해서는 다음의 값들이 필요하다.

• 측정 오차(easy to know)와 상태의 확률성(usually estimated)의 공분산 행렬인 R과 Q의 추정치
• 시간 0에서 이미 알고 있거나 추정된 상태 값 (x hat)_0 (y_0에 기반해 추정됨)
• 시간 t에 적용될 압력에 대한 사전 지식과, 그 압력이 상태에 어떤 영향을 미칠 것인가를 나타내는 행렬 B와 값 u_t
• 시간 단계에 따른 상태의 변화를 결정짓는 System dynamics의 knowledge, F
• 측정이 상태에 어떻게 의존하는 지에 대한 knowledge, A

a probabilistic perspective in expectation values, a least squares minimization problem, or as a maximum likelihood estimation problem 등 칼만 필터의 공식은 여러 방법으로 도출할 수 있다. (궁금하면 서치해보자)

 

---

Code for the Kalman Filter(with R)

 

• 사례: 오차가 쉽게 발생하는 센서로 뉴턴 역학을 따르는 물체를 추적한다.
  • 물체의 위치를 물체의 속도와 가속도에 대한 함수로 보는 Newton's laws of motion(뉴턴 운동 법칙)에 따라 시계열 생성
  • 가정: 물체의 움직임은 계속성을 갖지만, 측정이 이산적으로 분리됨
  • 가속도에 대한 계열 설정
  • 위치와 속도는 0에서 시작
  • (물리적으로 비현실적이지만) 각 시간 단계의 시작점에서 순간 가속도가 변하고, 다시 변하기 전까지는 지속적인 값을 가짐

 

## R

## 로켓은 시간 단계를 100번 거침. rocket will take 100 time steps
ts.length <- 100

## 가속도가 움직임을 주도. the acceleration will drive the motion
a <- rep(0.5, ts.length)

## 위치와 속도는 0에서 시작. position and velocity start at 0
x <- rep(0, ts.length)
v <- rep(0, ts.length)
for (ts in 2:ts.length) {
  x[ts] <- v[ts - 1] * 2 + x[ts - 1] + 1/2 * a[ts-1] ^ 2
  x[ts] <- x[ts] + rnorm(1, sd = 20) ## 확률적 요소. stochastic component
  v[ts] <- v[ts - 1] + 2 * a[ts-1]
}


par(mfrow = c(3, 1))
plot(x, main = "Position", type = 'l')
plot(v, main = "Velocity", type = 'l')
plot(a, main = "Acceleration", type = 'l')

 

Figure 1. 로켓의 위치, 속도, 가속도에 대한 그래프

Figure 1을 통해 가속도에 따른 움직임을 볼 수 있다.

여기서 모든 상태가 위치, 속도, 가속도의 세 변수로 묘사된다고 상정한다.

그러나, 우리에게 가용한 데이터는 오직 노이즈한 센서로부터 얻을 수 있는 물체의 위치뿐이다.

또한 선형적으로 증가하는 속도(velocity)에 영향을 주는 일정한 가속도(acceleration)가 위로 솟구치는 볼록 모양의 그래프(position)를 만들어낸다는 것을 알 수 있다.

 

 

다음 코드는 측정된 값이 실제 위치와 어떤 관계를 갖는지에 대한 그래프를 그린다(Figure 2).

여기서 변수 x가 바로 센서를 의미한다.

z <- x + rnorm(ts.length, sd = 300)
plot (x, ylim = range(c(x, z)))
lines(z)

 

Figure 2. 참인 위치(원)와 오차가 있는 측정(선) 간의 관계

(state transition equation에 넣은 노이즈로 인해 위치 x가 완벽한 볼록이 아니라는 사실에 유의)

 

 

이제 칼만 필터를 적용해보자. 먼저 앞에서 배운 식을 반영하는 일반적인 함수 코드를 작성한다.

## Kalman Filter
kalman.motion <- function(z, Q, R, A, H) {
  dimState = dim(Q)[1]
  
  xhatminus <- array(rep(0, ts.length * dimState),
                     c(ts.length, dimState)) # x hat-
  xhat <- array(rep(0, ts.length * dimState),
                c(ts.length, dimState)) # x hat
  Pminus <- array(rep(0, ts.length * dimState * dimState),
                  c(ts.length, dimState, dimState)) # P-
  P <- array(rep(0, ts.length * dimState * dimState),
             c(ts.length, dimState, dimState)) # P
  
  K <- array(rep(0, ts.length * dimState),
             c(ts.length, dimState)) # K, Kalman gain
  
  # 초기 추측 = 모든 지표는 0으로 시작함. intial guesses = starting at 0 for all metrics
  xhat[1, ] <- rep(0, dimState)
  P[1, , ] <- diag(dimState)
  
  for (k in 2:ts.length) {
    # Prediction
    xhatminus[k, ] <- A %*% matrix(xhat[k-1, ])
    Pminus[k, , ] <- A %*% P[k-1, , ] %*% t(A) + Q
    
    # Filtering
    K[k, ] <- Pminus[k, , ] %*% H %*%
      solve( t(H) %*% Pminus[k, , ] %*% H + R )
    xhat[k, ] <- xhatminus[k, ] + K[k, ] %*%
      (z[k]- t(H) %*% xhatminus[k, ])
    P[k, , ] <- (diag(dimState)-K[k,] %*% t(H)) %*% Pminus[k, , ]
  }
  
  ## 예측과 평활화된 값 모두 반환. we return both the forecast and the smoothed value
  return(list(xhat = xhat, xhatminus = xhatminus))
}

이제 이 함수를 적용하여 로켓의 위치만을 측정 가능하게 만들자(즉, 속도와 가속도는 포함하지 않음).

## noise parameters
R <- 10^2 ## 측정 분산 - 이 값은 측정 도구에 대해 알려진, 
## 물리적인 한계에 따라 설정되어야 한다.
## x에 더해진 노이즈와 일관성 있게 설정한다.

Q <- 10 ## 과정의 분산 - 일반적으로 성능의 최대화를 위해서
## 조정되어야 하는 하이퍼파라미터로 취급된다.

A <- matrix(1) ## x_t = A * x_(t-1) (이전 x가 나중의 x에 얼마나 영향을 미치는지. how prior x affects later x)
H <- matrix(1) ## y_t = H * x_t (상태를 측정으로 변환. translating state to measurement)

## 칼만 필터 method로 데이터를 넣고 돌림. run the data through the Kalman filtering method
xhat <- kalman.motion(z, diag(1) * Q, R, A, H)[[1]]

Figure 3. 측정 위치(Measured), 실제(참) 위치(Actual), 필터링된 추정된 위치(Filtered), 예측된 위치(Forecast)

• 측정 위치(Measured)
• 실제(참) 위치(Actual)
• 필터링된 추정된 위치(Filtered): 즉 시간 t에서의 위치와 시간 t의 측정을 결합한 최상의 추정
• 예측된 위치(Forecast): 즉 시간 t를 포함하지 않은 t-1까지의 측정에 알고 있는 system dynamics를 더한 것을, 시간 t의 위치에 결합한 최상의 추정

 

이처럼 칼만 필터는 측정 오차에서 많은 노이즈를 제거한다. 제거하는 정도는 측정 오차의 파라미터 R 값에 따라 달라진다. 이 값은 필터가 과거 대비 가장 최근의 데이터에 얼마나 큰 가중치를 부여하는지를 반영한다. 이렇게 필터는 꽤 만족스러운 데이터 예측을 만들어낸다. 특히, 예측과 실제 데이터 사이에 시간차가 없다는 사실에 주목하자. 즉, 직전 값에만 의존해서 현재 값을 예측했다는 것이다.

 

지금까지 간단한 칼만 필터 예제를 다루었다. 특히, 시스템의 내부 dynamics가 잘 정립되고 이해된 경우에 더욱 유용하며, 이상적인 도구가 될 수 있다(로켓 예제 등).

같은 내용이나 서로 다른 수치가 여러 장치에 의해 동시 다발적으로 측정되는 경우처럼, 여러 종류의 측정치를 다루는 경우 특히 유용하다.

 

칼만 필터의 한 가지 대단한 이점은 recursive(재귀성)이다. 즉 과정이 반복될 때마다 모든 데이터를 살펴볼 필요 없이, 각 시간 단계에서 그 이전 시간들로부터의 모든 정보가 가장 최근 상태, 공분산 추정이라는 몇 개의 파라미터로서 최선의 방법으로 결합되어 제공된다. 이러한 '요약 통계'를 사용하여, 지능적인 방식으로 update을 수행한다는 점이 칼만 필터의 장점이다. 이러한 장점은 계산의 자원이나 소요 시간이 매우 중요한 현실의 여러 적용 사례에 있어 칼만 필터의 유용성을 제고한다.

 

여기에서는 다루지 않았지만, 칼만 필터를 이용해 여러가지 유용한 확장을 할 수 있다. 칼만 필터를 평활화에 적용하는 것이 그 중 한 가지다. 즉 시간 t에서의 참 상태 추정에, 시간 t 전후 데이터를 모두 사용하여 최선의 추정을 한다는 것이다.

 

 

 

 

 

 

* code: https://github.com/inandout-kr/timeSeriesAnalysis/blob/main/KalmanFilter.R