13. Sequential Data

순차적 데이터를 다루는 방법을 알아보자.

13.1. Markov Models

마르코프 연쇄의 기본 발상은 X_{t}가 미래를 예측하는 데 있어 충족 통계량이란 것이다. 이산적 시간 단계를 가정한다면 결합분포를 다음과 같이 쓸 수 있다:

p(X_{1:T}) = p(X_{1}) \prod_{t=2}^{T} p(X_{t} | X_{t-1})

이를 마르코프 연쇄 또는 마르코프 모델이라 한다.

13.2. Hidden Markov Models

은닉 마르코프 모델(HMM)은 은닉 상태 z_{t} \in \{1, \cdots, K \}과 그에 대한 관찰 p(\mathbf{x}_{t} | z_{t})로 이루어진 이산 시간의 이산 상태 마르코프 연쇄이다. 그에 해당하는 결합분포는 p(\mathbf{z}_{1 : T} | \mathbf{x}_{1:T}) = p(\mathbf{z}_{1:T}) p(\mathbf{x}_{1:T} | \mathbf{z}_{1: T}) = [p(z_{1}) \prod_{t=2}^{T} p(z_{t} | z_{t-1})] [\prod_{t=1}^{T} p(\mathbf{x}_{t} | z_{t})]의 형태가 된다.

은닉 마르코프 모델의 은닉 상태는 이산적이지만 그에 대한 관찰은 이산적일 수도 있고 연속적일 수도 있다. 이산적인 경우에는 관찰 모델을 관찰 행렬 p(\mathbf{x}_{t} = l | z_{t} = k, \mathbf{\theta}) = B_{kl}로 잡는다. 연속적인 경우에는 관찰 모델을 조건부 가우시안 p(\mathbf{x}_{t} | z_{t} = k, \mathbf{\theta}) = \mathcal{N}(\mathbf{x}_{t} | \mathbf{\mu}_{k}, \mathbf{\Sigma}_{k})로 잡는다. 이는 가우시안 혼합 모델과 비슷하지만 클러스터가 마르코프 역학계를 따른다는 것이 다르다. 그래서 이는 마르코프 치환 모델이라고도 불린다.

13.2.1. Maximum likelihood for the HMM

기대값 최대화 알고리즘(Baum-Welch)를 사용한다.

기대완전데이터 로그가능도는 다음과 같다:

Q(\mathbf{\theta}, \mathbf{\theta}_{\mathrm{old}}) = \sum_{k=1}^{K} \mathbb{E}[N_{k}^{1}] \log \pi_{k} + \sum_{j=1}^{K} \sum_{k=1}^{K} \mathbb{E}[N_{jk}] \log A_{jk} + \sum_{i=1}^{N} \sum_{t=1}^{T_{i}} \sum_{k=1}^{K} p(z_{t} = k | \mathbf{x}_{i}, \mathbf{\theta}^{\mathrm{old}}) \log p(\mathbf{x}_{i, t} | \mathbf{\phi}_{k})

\mathbb{E}[N_{k}^{1}] = \sum_{i=1}^{N} p(z_{i1} = k | \mathbf{x}_{i}, \mathbf{\theta}^{\mathrm{old}})

\mathbb{E}[N_{jk}] = \sum_{i=1}^{N} \sum_{t=2}^{T_{i}} p(z_{i, t - 1} = j, z_{i, t} = k | \mathbf{x}_{i}, \mathbf{\theta}^{\mathrm{old}})

\mathbb{E}[N_{j}] = \sum_{i=1}^{N} \sum_{t=1}^{T_{i}} p(z_{i, t} = j | \mathbf{x}_{i}, \mathbf{\theta}^{\mathrm{old}})

이 충족 통계량은 전방-후방 알고리즘을 통해 \gamma_{i, t}(j) = p(z_{t} = j | \mathbf{x}_{i, 1 : T_{i}}, \mathbf{\theta}), \xi_{i, t}(j, k) = p(z_{t-1} = j, z_{t} = k | \mathbf{x}_{i, 1 : T_{i}}, \mathbf{\theta})를 계산함으로써 계산할 수 있다.

\mathbf{A}\mathbf{\pi}의 M 단계는 기대 등장빈도 \hat{A}_{jk} = \frac{\mathbb{E}[N_{jk}]}{\sum_{k^{\prime}} \mathbb{E} [N_{jk^{\prime}}]}, \hat{\pi}_{k} = \frac{\mathbb{E}[N_{k}^{1}]}{N}이다. 이는 직관적으로 상태 j에서 k로 전이되는 횟수의 기대값을 더한 뒤 상태 j에서 다른 상태로 전이되는 횟수로 나눈 것과 같다.

13.2.2. The forward-backward algorithm

위에서는 온라인 추론을 통해 필터링된 주변분포 p(z_{t} = j | \mathbf{x}_{1:t})를 계산하였다. 여기서는 오프라인 추론을 통해 다듬질된 주변분포 p(z_{t} = j | \mathbf{x}_{1:T})을 계산한다. 기본 아이디어는 연쇄를 과거/미래 두 부분으로 나누는 것이다.

p(z_{t} = j | \mathbf{x}_{1:T}) \propto p(z_{t} = j, \mathbf{x}_{t+1:T} | \mathbf{x}_{1:T}) \propto p(z_{t} = j | \mathbf{x}_{1:t}) p(\mathbf{x}_{t+1 : T} | z_{t} = j)

이 때 \alpha_{t}(j) = p(z_{t} = j | \mathbf{x}_{1:t})를 필터링된 믿음 상태, \beta_{t}(j) = p(\mathbf{x}_{t+1:T} | z_{t} = j)을 은닉 상태에 대해 조건을 건 미래 증거 조건 가능도, \gamma_{t}(j) = p(z_{t} = j | \mathbf{x}_{1:T})을 다듬질된 사후주변분포라 하자. 그러면 \gamma_{t}(j) \propto \alpha_{t}(j) \beta_{t}(j)이 성립한다.

위에서 \mathbf{\alpha}_{t}를 재귀적으로 구했듯이 \mathbf{\beta}_{t}도 재귀적으로 구할 수 있는데 \mathbf{\beta}_{t-1} = \mathbf{\Psi}(\mathbf{\psi}_{t} \odot \mathbf{\beta}_{t})이 된다. 기본 케이스는 \beta_{T}(i) = 1로, 이벤트가 발생하지 않았을 때의 확률이다. 전방-후방 메시지를 계산한 뒤에는 \gamma_{t}(j) \propto \alpha_{t}(j) \beta_{t}(j)로 사후주변분포를 계산할 수 있는데, 이를 전방-후방 알고리즘이라 한다.

13.2.3. The sum-product algorithm for the HMM

더 자세히 알아보자면, 정점 t의 믿음 상태를 계산하고 싶다고 할 때, 최초에는 트리에서 정점 t의 후손들의 믿음 상태에 대해서만 조건을 건다. 즉, \mathrm{bel}_{t}^{-}(x_{t}) = p(x_{t} | \mathbf{v}_{t}^{-})로 놓는다. 이제 t의 자식 노드 s에 대해 s를 루트로 하는 부분 트리에 대해 t가 알아야 하는 정보를 메시지 m_{s \to t}^{-}(x_{t}) = p(x_{t} | \mathbf{v}_{st}^{-})로 두면, 최종 믿음 상태는 다음과 같다.

\mathrm{bel}_{t}^{-}(x_{t}) = p(x_{t} | \mathbf{v}_{t}^{-}) = \frac{1}{Z_{t}} \phi_{t}(x_{t}) \prod_{c \in \mathrm{ch}(t)} m_{c \to t}^{-} (x_{t})

즉, 노드의 국소적 증거와 자식 노드로부터 오는 정보들을 전부 곱한 뒤 표준화하는 것이다. 그러면 자식 노드로부터 오는 정보 메시지는 어떻게 계산해야 할까? 이는 다음과 같다.

m_{s \ to t}^{-}(x_{t}) = \sum_{x_{s}} \phi_{st}(x_{s}, x_{t}) \mathrm{bel}_{s}^{-}(x_{s})

즉, 정점 s에 대한 믿음 상태를 정점 t에 대한 믿음 상태와 s와 t를 잇는 간선의 잠재 함수에 대한 식으로 치환한 것이다. 이를 루트에 도달할 때까지 반복한다. 루트에서는 트리 내의 모든 증거량을 관측하였으므로 루트의 국소적 믿음 상태는 다음을 이용해 계산할 수 있다.

\mathrm{bel}_{r}(x_{r}) = p(x_{r} | \mathbf{v}) = p(x_{t} | \mathbf{v}_{r}^{-}) \propto \phi_{r}(x_{r}) \prod_{c \in \mathrm{ch}(r)} m_{c \to r}^{-}(x_{r}).

이는 알고리즘의 상승 방향 부분을 완성시킨다. 부가적으로, 증거량에 대한 확률을 표준화 상수의 곱으로 계산할 수도 있다: p(\mathbf{v}) = \prod_{t} Z_{t}

이제 루트로부터 메시지를 내려보내 보자. 정점 t에서 자식 s로 보내는 메시지를 m_{t \to s}^{+} (x_{s}) = p(x_{t} | \mathbf{v}_{st}^{+})라 하면 \mathrm{bel}_{s}(x_{s}) = p(x_{s} | \mathbf{v}) \propto \mathrm{bel}_{s}^{-}(x_{s}) \prod_{t \in \mathrm{pa}(s)} m_{t \to s}^{+}(x_{s})이 된다.

그러면 아래쪽 방향으로 보내는 메시지는 어떻게 계산할까? 이는 다음과 같다.

m_{t \to s}^{+}(x_{s}) = p(x_{s} | \mathbf{v}_{st}^{+}) = \sum_{x_{t}} \phi_{st}(x_{s}, x_{t}) \frac{\mathrm{bel}_{t}(x_{t})}{m_{s \to t}^{-}(x_{t})} = \sum_{x_{t}} \phi_{st}(x_{s}, x_{t}) \prod_{c \in \mathrm{ch}(t), c \neq s} m_{c \to t}^{-} (x_{t}) \prod_{p \in \mathrm{pa}(t)} m_{p \to t}^{+}(x_{t})

즉, 자식 s를 제외한 모든 노드로부터 정점 t로 오는 모든 메시지들을 곱한 뒤 이를 간선 s – t에 대한 잠재 함수로 전파한다. 연쇄 구조인 경우에는 (즉, 트리에서 자식과 부모가 하나인 경우) 위는 다음과 같다.

m_{t \to s}^{+}(x_{s}) = p(x_{s} | \mathbf{v}_{st}^{+}) = \sum_{x_{t}} \phi_{st}(x_{s}, x_{t}) \phi_{t}(x_{t}) m_{p \to t}^{+}(x_{t})

아래쪽 방향으로 보내는 메시지를 계산할 때 모든 메시지를 하나 제외하고 곱하는 방법을 합-곱 방법이라 한다. 이는 전방-후방 알고리즘의 후방 부분과 비슷하다.

13.2.4. Scaling factors

전방-후방 알고리즘을 실제로 쓸 때에는 언더플로우를 방지하기 위해 로그를 취하거나, 항들을 재스케일링해야 한다. 이 스케일 인자는 기대값 최대화 과정에서 약분된다.

13.2.5. The Viterbi algorithm

비터비 알고리즘은 연쇄 구조 그래프 모델에서 가장 확률이 높은 서열을 계산한다. 즉, \mathbf{z}^{\ast} = \mathrm{argmax}_{\mathbf{z}_{1:T}} p(\mathbf{z}_{1:T} | \mathbf{x}_{1:T})을 계산한다. \delta_{j} = \max_{z_{1}, \cdots, z_{t-1}} p(\mathbf{z}_{1 : t-1}, z_{t} = j | \mathbf{x}_{1:t})를 가장 가능성 높은 경로를 탔을 때 시간 t에서 상태 j에 도달할 확률이라 하자. 이는 시간 t-1에서 상태 i에 도달할 확률에 상태 i에서 상태 j로 전이할 확률을 곱한 것이므로 \delta_{t}(j) = \max{i} \delta_{t-1}(i) \phi(i, j) \phi_{t}(j)이 된다. 가장 확률이 높은 직전 상태 a_{t}(j) = \mathrm{argmax}_{i} \delta_{t-1}(i) \psi(i, j) \phi_{t}(j)도 저장해 놓아야 한다. 초기값은 \delta_{1}(j) = \pi_{j} \phi_{1}(j)로 선택하고 최종 상태는 z_{T}^{\ast} = \mathrm{argmax}_{i} \delta_{T}(i)이 된다. 이 후에 역추적으로 z_{t}^{\ast} = a_{t+1}(z_{t+1}^{\ast})을 통해 계산해 나가면 된다.

13.2.6. Extensions of the hidden Markov model

은닉 마르코프 모델을 보완하기 위한 여러 변형들이 존재한다. 자가회귀적 은닉 마르코프 모델, 인자 은닉 마르코프 모델 등.

13.3. Linear Dynamical Systems

상태 공간 모델의 중요한 특수 케이스는 모든 조건부확률분포가 선형 가우시안인 경우이다. 다시 말해,

\mathbf{z}_{t} = \mathbf{A}_{t} \mathbf{z}_{t-1} + \mathbf{B}_{t} \mathbf{u}_{t} + \mathbf{\epsilon}_{t}

\mathbf{y}_{t} = \mathbf{C}_{t} \mathbf{z}_{t} + \mathbf{D}_{t} \mathbf{u}_{t} + \mathbf{\delta}_{t}

\mathbf{\epsilon}_{t} \sim \mathcal{N}(\mathbf{0}, \mathbf{Q}_{t})

\mathbf{\delta}_{t} \sim \mathcal{N}(\mathbf{0}, \mathbf{R}_{t})

인 경우를 말한다. 이를 선형-가우시안 상태 공간 모델 (LG-SSM) 또는 선형 동역학계 (LDS) 라고도 한다. 만약 모델 매개변수 \mathbf{\theta}_{t} = (\mathbf{A}_{t}, \mathbf{B}_{t}, \mathbf{C}_{t}, \mathbf{D}_{t}, \mathbf{Q}_{t}, \mathbf{R}_{t})이 시간에 독립적이라면, 이 모델은 정적이라 한다. 이 케이스가 중요한 이유는 정확한 추론이 가능하기 때문인데, 최초의 믿음 상태가 가우시안 p(\mathbf{z}_{1}) = \mathcal{N}(\mathbf{\mu}_{1 | 0}, \mathbf{\Sigma}_{1 |0})이라면 모든 믿음 상태 또한 가우시안 p(\mathbf{z}_{t} | \mathbf{y}_{1 : t}) = \mathcal{N}(\mathbf{\mu}_{t | t}, \mathbf{\Sigma}_{t | t})이 된다. 이는 칼만 필터로 계산할 수 있다.

13.3.1. Inference in LDS

칼만 거르기는 선형 가우시안 상태 공간 모델에 대한 정확한 베이지안 거르기 알고리즘이다. 시점 t에서의 주변사후분포를 p(\mathbf{z}_{t} | \mathbf{y}_{1 : t}, \mathbf{u}_{1 : t}) = \mathcal{N}(\mathbf{z}_{t} | \mathbf{\mu}_{t}, \mathbf{\Sigma}_{t})로 나타내면 모든 것이 가우시안이므로 예측과 업데이트를 모두 닫힌 형태로 표현할 수 있다. 이는 은닉 마르코프 모델과 비슷하다.

예측 단계는 다음과 같다.

p(\mathbf{z}_{t} | \mathbf{y}_{1 : t-1}, \mathbf{u}_{1 : t}) = \int \mathcal{N}(\mathbf{z}_{t} | \mathbf{A}_{t} \mathbf{z}_{t-1} + \mathbf{B}_{t} \mathbf{u}_{t}, \mathbf{Q}_{t}) \mathcal{N}(\mathbf{z}_{t-1} | \mathbf{\mu}_{t-1}, \mathbf{\Sigma}_{t-1}) d \mathbf{z}_{t-1} = \mathcal{N}(\mathbf{z}_{t} | \mathbf{\mu}_{t | t - 1}, \mathbf{\Sigma}_{t | t - 1})

\mathbf{\mu}_{t | t - 1} = \mathbf{A}_{t} \mathbf{\mu}_{t-1} + \mathbf{B}_{t} \mathbf{u}_{t}

\mathbf{\Sigma}_{t | t - 1} = \mathbf{A}_{t} \mathbf{\Sigma}_{t - 1} \mathbf{A}_{t}^{T} + \mathbf{Q}_{t}

로그주변가능도도 계산할 수 있다.

\log p(\mathbf{y}_{1 : T} | \mathbf{u}_{1 : T}) = \sum_{t} \log p(\mathbf{y}_{t} | \mathbf{y}_{1 : t - 1}, \mathbf{u}_{1 : t}) = \sum_{t} \log \mathcal{N}(\mathbf{y}_{t} | \mathbf{C}_{t} \mathbf{\mu}_{t | t- 1}, \mathbf{S}_{t})

한 단계 전 사후예측분포는 다음과 같으며, 시계열 예측에 유용하다.

p(\mathbf{y}_{t} | \mathbf{y}_{1 : t - 1}, \mathbf{u}_{1 : t}) = \mathcal{N}(\mathbf{y}_{t} | \mathbf{C} \mathbf{\mu}_{t | t - 1}, \mathbf{C} \mathbf{\Sigma}_{t | t - 1}\mathbf{C}^{T} + \mathbf{R})

13.3.2. Learning in LDS

출력 서열만 관측할 수 있다면, 매개변수에 대한 최대가능도추정이나 최대사후확률추정은 기대값 최대화 알고리즘으로 이루어진다. 이는 은닉 마르코프 모델에서의 바움-웰치 알고리즘과 비슷하지만, E 단계에서 전방-후방 알고리즘 대신 칼만 다듬질 알고리즘을 사용하므로 M 단계에서의 계산도 달라진다는 점이 다르다.

13.3.3. Extensions of LDS

노이즈가 가우시안이고 전이 모델 g와 관측 모델 h는 비선형이되 미분가능한 다음의 비선형 모델을 가정하자.

\mathbf{z}_{t} = g(\mathbf{u}_{t}, \mathbf{z}_{t - 1}) + \mathcal{N}(\mathbf{0}, \mathbf{Q}_{t})

\mathbf{y}_{t} = h(\mathbf{z}_{t}) + \mathcal{N}(\mathbf{0}, \mathbf{R}_{t})

여기서는 사후분포를 단일 가우시안으로 근사할 수 있는 경우에 집중해보도록 하자. 확장된 칼만 거르기(EKF)는 이러한 형태의 비선형 가우시안 동역학계에 쓰일 수 있다. 기본 발상은 g와 h를 일차 테일러 급수 전개 근사를 이용해 이전 상태 추정값에 대한 선형 함수로 근사하는 것이다. 이후 표준 칼만 거르기 식들을 적용하면 된다. (노이즈 \mathbf{Q}, \mathbf{R}은 변하지 않는다고 가정한다.) 즉, 정적인 비선형 동역학계를 비-정적인 선형 동역학계로 근사하는 것이다.

\mathbf{H}_{t}\mathbf{h}의 자코비안을 사전최빈값에서 구한 값이라 할 때, 측정 모델은 다음과 같이 근사된다:

p(\mathbf{y}_{t} | \mathbf{z}_{t}) \simeq \mathcal{N}(\mathbf{y}_{t} | \mathbf{h}(\mathbf{\mu}_{t | t - 1}) + \mathbf{H}_{t} (\mathbf{y}_{t} - \mathbf{\mu}_{t | t - 1}), \mathbf{R}_{t})

\mathbf{H}_{t} = \mathbf{H} |_{\mathbf{z} = \mathbf{\mu}_{t | t - 1}}

\mathbf{G}_{t}\mathbf{g}의 자코비안을 사전최빈값에서 구한 값이라 할 때, 계 모델은 다음과 같이 근사된다:

p(\mathbf{z}_{t} | \mathbf{z}_{t - 1}, \mathbf{u}_{t}) \simeq \mathcal{N}(\mathbf{z}_{t} | \mathbf{g}(\mathbf{u}_{t}, \mathbf{\mu}_{t - 1}) + \mathbf{G}_{t}(\mathbf{z}_{t - 1} - \mathbf{\mu}_{t - 1}), \mathbf{Q}_{t})

\mathbf{G}_{t} = \mathbf{G}(\mathbf{u}_{t}) |_{\mathbf{z} = \mathbf{\mu}_{t | t - 1}}

이를 통해 칼만 거르기를 적용하면 사후분포를 다음과 같이 구할 수 있다.

\mathbf{\mu}_{t | t - 1} = \mathbf{g}(\mathbf{u}_{t}, \mathbf{\mu}_{t - 1})

\mathbf{V}_{t | t - 1} = \mathbf{G}_{t} \mathbf{V}_{t - 1} \mathbf{G}_{t}^{T} + \mathbf{Q}_{t}

\mathbf{K}_{t} = \mathbf{V}_{t | t - 1}\mathbf{H}_{t}^{T} (\mathbf{H}_{t} \mathbf{V}_{t | t - 1}\mathbf{H}_{t}^{T} + \mathbf{R}_{t})^{-1}

\mathbf{\mu}_{t} = \mathbf{\mu}_{t | t - 1} + \mathbf{K}_{t}(\mathbf{y}_{t} - \mathbf{h}(\mathbf{\mu}_{t | t - 1}))

\mathbf{V}_{t} = (\mathbf{I} - \mathbf{K}_{t} \mathbf{H}_{t}) \mathbf{V}_{t | t - 1}

일반 칼만 거르기와 다른 점은 상태를 예측할 때 \mathbf{A}_{t} \mathbf{\mu}_{t - 1} + \mathbf{B}_{t} \mathbf{\mu}_{t} 대신 \mathbf{g}(\mathbf{u}_{t}, \mathbf{\mu}_{t - 1})을 사용한다는 것이고 측정을 업데이트할 때 \mathbf{C}_{t} \mathbf{\mu}_{t | t - 1} 대신 h(\mathbf{\mu}_{t | t - 1})을 사용한다는 것이다. \mathbf{\mu}_{t | t - 1} 대신 \mathbf{\mu}_{t}를 중심으로 선형 근사를 수행해 성능을 향상시킬 수 있는데 이를 반복 확장된 칼만 거르기라고 한다. 더 느리지만 더 좋은 성능을 낸다.

확장된 칼만 거르기가 잘 동작하지 않는 두 가지 경우가 있는데 첫째는 사전 공분산이 클 경우 사전분포가 넓어지기 때문에 선형 근사가 평균에서 많이 떨어진 곳으로 잡히게 된다. 둘째는 현재 평균 근처에서 함수가 심하게 비선형적일 경우 당연히 선형 근사가 잘 동작하지 않는다. 이는 무취 칼만 거르기라는 알고리즘으로 보정할 수 있다.

13.3.4. Particle filters

입자 거르기(PF)는 몬테 카를로, 또는 시뮬레이션 기반 알고리즘으로 재귀적 베이지안 추론을 수행한다. 즉, 예측-업데이트 루프를 근사한다. 기본 발상은 믿음 상태를 입자들의 가중치 집합으로 근사하는 것이다.

p(\mathbf{z}_{1 : t} | \mathbf{y}_{1 : t}) \simeq \sum_{s=1}^{S} \hat{w}_{t, s} \delta_{\mathbf{z}_{1 : t, s}}(\mathbf{z}_{1 : t})

이로부터 가장 최근 상태의 주변분포는 p(\mathbf{z}_{t} | \mathbf{y}_{1 : t}) \simeq \sum_{s=1}^{S} \tilde{w}_{t, s} \delta_{\mathbf{z}_{t, s}}(\mathbf{z}_{t})로 근사할 수 있다.