Bayes and Particle Filters·Lecture 7 of 7

L7 · Resampling and SIR

Updated 2026-09-15·KO

시리즈: Bayes and Particle Filters · 이전: L6 · concept: Particle Filter

이 강의는 weighted particle을 시간에 따라 계속 들고 가면 weight가 소수의 particle에 몰리는 이유와, resampling이 그중 무엇을 고치고 무엇을 고치지 못하는지 답합니다. 마지막에 propagate, weight, normalize, resample을 하나의 SIR 알고리즘으로 묶고 Bayes filter의 식과 한 줄씩 대응시킵니다.

형광펜: 정의·핵심 · 직관·메커니즘 · 조건·주의

목차


1. Weight를 계속 들고 가면 생기는 일

출발점

L6 §6에서 transition model을 proposal로 쓰면 weight가 likelihood 하나로 정해진다는 것을 보았습니다.

wt(i)∝p(yt∣xt(i))w_t^{(i)} \propto p(y_t \mid x_t^{(i)})

그 유도에는 전제가 있었습니다. 시각 t−1t-1의 particle이 모두 같은 weight를 가지거나, 조상을 weight에 비례해서 골라야 새 particle이 mixture ∑jwˉt−1(j) p(xt∣xt−1(j))\sum_j \bar w_{t-1}^{(j)}\,p(x_t \mid x_{t-1}^{(j)})에서 뽑힌 것이 됩니다.

이제 그 전제를 빼 봅니다. 시각 t−1t-1에 weight가 서로 다른 particle 묶음 {xt−1(i), wt−1(i)}\{x_{t-1}^{(i)},\,w_{t-1}^{(i)}\}이 있고, 각 particle은 자기 조상에서만 transition을 한 번 뽑습니다. 이렇게 resampling 없이 propagate와 weight만 반복하는 방법을 sequential importance sampling (SIS)이라고 합니다. 이때 weight는 다음처럼 갱신됩니다.

wt(i)∝wt−1(i) p(yt∣xt(i))w_t^{(i)} \propto w_{t-1}^{(i)}\,p(y_t \mid x_t^{(i)})

직관은 "지금까지의 measurement를 얼마나 잘 설명했는가"에 "이번 measurement를 얼마나 잘 설명하는가"를 곱한다는 것입니다. 하지만 L6 §6의 논증은 이 경우에 그대로 쓸 수 없으므로 따로 유도합니다.

Target을 path posterior로 바꾼다

Particle ii가 자기 조상만 따라가면, 시각 tt까지 만든 것은 점 하나가 아니라 trajectory x0:t(i)=(x0(i),…,xt(i))x_{0:t}^{(i)}=(x_0^{(i)},\dots,x_t^{(i)})입니다. Transition model을 proposal로 쓰면 이 trajectory가 뽑힌 density는 한 step씩 곱으로 쌓입니다.

qt(x0:t)=qt−1(x0:t−1)  p(xt∣xt−1),q0(x0)=p(x0)q_t(x_{0:t}) = q_{t-1}(x_{0:t-1})\;p(x_t \mid x_{t-1}),\qquad q_0(x_0)=p(x_0)

"이전 path는 그대로 두고 끝에 transition sample 하나를 붙인다"는 뜻입니다. 그래서 importance sampling도 trajectory 위에서 합니다. SIS의 target은 filtering distribution p(xt∣y1:t)p(x_t \mid y_{1:t})가 아니라 path posterior p(x0:t∣y1:t)p(x_{0:t} \mid y_{1:t})입니다. p(x0:t∣y1:t)p(x_{0:t} \mid y_{1:t})는 "y1y_1부터 yty_t까지 본 뒤 trajectory x0,…,xtx_0,\dots,x_t 전체에 대한 posterior"라고 읽습니다.

Target을 바꾸는 이유가 하나 더 있습니다. L6 §4에서 importance sampling은 unnormalized target을 점마다 계산할 수 있어야 한다고 했습니다. Exact filtering distribution은 p(yt∣xt) bel‾t(xt)p(y_t\mid x_t)\,\overline{\mathrm{bel}}_t(x_t)에 비례하고, bel‾t\overline{\mathrm{bel}}_t 자체가 적분이라 점마다 계산되지 않습니다.

L6 §6처럼 bel‾t\overline{\mathrm{bel}}_t를 particle mixture ∑jwˉt−1(j) p(xt∣xt−1(j))\sum_j \bar w_{t-1}^{(j)}\,p(x_t\mid x_{t-1}^{(j)})로 바꾸면 점마다 계산할 수는 있습니다. 하지만 weight가 서로 다른 채로 각 particle이 자기 조상만 따라가면, 조상이 weight와 상관없이 한 번씩 쓰이므로 proposal이 그 mixture가 아닙니다. 그래서 L6 §6의 약분이 성립하지 않습니다. 반면 path posterior는 L3 §5의 분해에 따라 p(x0)∏k=1tp(xk∣xk−1) p(yk∣xk)p(x_0)\prod_{k=1}^t p(x_k\mid x_{k-1})\,p(y_k\mid x_k)에 비례하므로, model의 factor를 곱하기만 하면 점마다 계산됩니다.

Path posterior의 재귀식

y1:t=(y1:t−1, yt)y_{1:t}=(y_{1:t-1},\,y_t)로 나누고, y1:t−1y_{1:t-1}이 주어진 세계에서 식을 전개합니다.

p(x0:t∣y1:t)=p(x0:t, yt∣y1:t−1)p(yt∣y1:t−1)(conditional probability given y1:t−1)=p(yt∣x0:t, y1:t−1)  p(x0:t∣y1:t−1)p(yt∣y1:t−1)(product rule)=p(yt∣xt)  p(x0:t∣y1:t−1)p(yt∣y1:t−1)(measurement CI)=p(yt∣xt)  p(xt∣x0:t−1, y1:t−1)  p(x0:t−1∣y1:t−1)p(yt∣y1:t−1)(product rule)=p(yt∣xt)  p(xt∣xt−1)  p(x0:t−1∣y1:t−1)p(yt∣y1:t−1)(Markov assumption)\begin{aligned} p(x_{0:t} \mid y_{1:t}) &= \frac{p(x_{0:t},\,y_t \mid y_{1:t-1})}{p(y_t \mid y_{1:t-1})} && \text{(conditional probability given } y_{1:t-1})\\ &= \frac{p(y_t \mid x_{0:t},\,y_{1:t-1})\;p(x_{0:t} \mid y_{1:t-1})}{p(y_t \mid y_{1:t-1})} && \text{(product rule)}\\ &= \frac{p(y_t \mid x_t)\;p(x_{0:t} \mid y_{1:t-1})}{p(y_t \mid y_{1:t-1})} && \text{(measurement CI)}\\ &= \frac{p(y_t \mid x_t)\;p(x_t \mid x_{0:t-1},\,y_{1:t-1})\;p(x_{0:t-1} \mid y_{1:t-1})}{p(y_t \mid y_{1:t-1})} && \text{(product rule)}\\ &= \frac{p(y_t \mid x_t)\;p(x_t \mid x_{t-1})\;p(x_{0:t-1} \mid y_{1:t-1})}{p(y_t \mid y_{1:t-1})} && \text{(Markov assumption)} \end{aligned}

셋째 줄은 L3의 measurement 가정 p(yt∣x0:t,y1:t−1)=p(yt∣xt)p(y_t\mid x_{0:t},y_{1:t-1})=p(y_t\mid x_t)이고, 다섯째 줄은 full form의 Markov assumption p(xt∣x0:t−1,y1:t−1)=p(xt∣xt−1)p(x_t\mid x_{0:t-1},y_{1:t-1})=p(x_t\mid x_{t-1})입니다. 두 가정이 각각 한 번씩 쓰였습니다. 분모 p(yt∣y1:t−1)p(y_t\mid y_{1:t-1})은 trajectory에 따라 달라지지 않으므로 비례식으로 쓰면 사라집니다.

Weight의 재귀식

Path 위의 importance weight는 target과 proposal의 비입니다. 모든 particle에 공통인 상수는 ∝\propto에 흡수합니다.

wt(i)∝p(x0:t(i)∣y1:t)qt(x0:t(i))(importance weight on path space)∝p(yt∣xt(i))  p(xt(i)∣xt−1(i))  p(x0:t−1(i)∣y1:t−1)qt−1(x0:t−1(i))  p(xt(i)∣xt−1(i))(recursion above, proposal)=p(yt∣xt(i))  p(x0:t−1(i)∣y1:t−1)qt−1(x0:t−1(i))(transition cancels)∝p(yt∣xt(i))  wt−1(i)(definition of wt−1(i))\begin{aligned} w_t^{(i)} &\propto \frac{p(x_{0:t}^{(i)} \mid y_{1:t})}{q_t(x_{0:t}^{(i)})} && \text{(importance weight on path space)}\\ &\propto \frac{p(y_t \mid x_t^{(i)})\;p(x_t^{(i)} \mid x_{t-1}^{(i)})\;p(x_{0:t-1}^{(i)} \mid y_{1:t-1})}{q_{t-1}(x_{0:t-1}^{(i)})\;p(x_t^{(i)} \mid x_{t-1}^{(i)})} && \text{(recursion above, proposal)}\\ &= p(y_t \mid x_t^{(i)})\;\frac{p(x_{0:t-1}^{(i)} \mid y_{1:t-1})}{q_{t-1}(x_{0:t-1}^{(i)})} && \text{(transition cancels)}\\ &\propto p(y_t \mid x_t^{(i)})\;w_{t-1}^{(i)} && \text{(definition of } w_{t-1}^{(i)}) \end{aligned}
  wt(i)∝wt−1(i)  p(yt∣xt(i))  \boxed{\;w_t^{(i)} \propto w_{t-1}^{(i)}\;p(y_t \mid x_t^{(i)})\;}

Transition density는 target에도 proposal에도 같은 값으로 들어가서 약분됩니다. 그래서 transition density의 값을 계산할 필요가 없고, transition에서 sample을 뽑을 수만 있으면 됩니다. w0(i)=1/Nw_0^{(i)}=1/N에서 시작해 이 식을 tt번 풀어 쓰면 다음과 같습니다.

wt(i)∝∏k=1tp(yk∣xk(i))w_t^{(i)} \propto \prod_{k=1}^{t} p(y_k \mid x_k^{(i)})

SIS의 weight는 그 trajectory가 받은 likelihood 전체의 곱입니다. 이 곱이 §2의 주제입니다.

Filtering distribution 읽어내기

Normalized weight를 wˉt(i)=wt(i)/∑jwt(j)\bar w_t^{(i)}=w_t^{(i)}/\sum_j w_t^{(j)}로 두면, weighted path 묶음은 expectation을 근사한다는 뜻에서 path posterior를 나타냅니다. Path 위의 Dirac delta는 시각마다의 delta의 곱 δ(x0:t−x0:t(i))=∏k=0tδ(xk−xk(i))\delta(x_{0:t}-x_{0:t}^{(i)})=\prod_{k=0}^t\delta(x_k-x_k^{(i)})입니다. 앞부분을 marginalization으로 없앱니다.

p^(xt∣y1:t)=∫∑i=1Nwˉt(i) δ(x0:t−x0:t(i)) dx0:t−1(marginalization)=∑i=1Nwˉt(i) δ(xt−xt(i))(∫δ(xk−xk(i)) dxk=1 for k<t)\begin{aligned} \hat p(x_t \mid y_{1:t}) &= \int \sum_{i=1}^N \bar w_t^{(i)}\,\delta\big(x_{0:t}-x_{0:t}^{(i)}\big)\,dx_{0:t-1} && \text{(marginalization)}\\ &= \sum_{i=1}^N \bar w_t^{(i)}\,\delta\big(x_t-x_t^{(i)}\big) && \Big(\textstyle\int \delta(x_k-x_k^{(i)})\,dx_k=1 \text{ for } k<t\Big) \end{aligned}

Path 위에서 구한 weight를 그대로 두고 마지막 값 xt(i)x_t^{(i)}만 보면 filtering distribution의 근사가 됩니다. Weight 재귀식에도 xt(i)x_t^{(i)}와 wt−1(i)w_{t-1}^{(i)}만 나오므로, 실제 구현은 path 전체를 저장하지 않습니다.

매 step resampling하면

§4의 resampling을 매 step 하면, resampling된 path들은 resampling 전 묶음이 만든 weighted empirical path distribution p^(x0:t−1∣y1:t−1)=∑iwˉt−1(i) δ(x0:t−1−x0:t−1(i))\hat p(x_{0:t-1}\mid y_{1:t-1})=\sum_i \bar w_{t-1}^{(i)}\,\delta\big(x_{0:t-1}-x_{0:t-1}^{(i)}\big)에서 뽑은 같은 weight의 sample입니다. 시각 t−1t-1의 target으로 p(x0:t−1∣y1:t−1)p(x_{0:t-1}\mid y_{1:t-1}) 대신 이 근사를 쓰면, sample이 그 target에서 직접 뽑혔으므로 wt−1(i)=1/Nw_{t-1}^{(i)}=1/N입니다. 이 값은 모든 particle에 공통이므로 ∝\propto에 흡수되고 wt(i)∝p(yt∣xt(i))w_t^{(i)}\propto p(y_t\mid x_t^{(i)})만 남습니다. L6 §6의 mixture 논증과 여기의 path space 논증은 서로 다른 논증이지만, 매 step resampling하면 둘 다 근사 target에 대해 같은 식이 나옵니다.


2. Degeneracy는 왜 생기는가

Degeneracy라는 상태

Particle 5개의 normalized weight가 처음에는 [0.2, 0.2, 0.2, 0.2, 0.2][0.2,\,0.2,\,0.2,\,0.2,\,0.2]였는데, SIS로 여러 step을 지나고 나니 [0.0001, 0.0002, 0.9989, 0.0005, 0.0003][0.0001,\,0.0002,\,0.9989,\,0.0005,\,0.0003]이 되었다고 합시다. Particle은 여전히 5개이고 5개 모두 propagate하고 likelihood를 계산합니다. 하지만 weighted estimate ∑iwˉ(i)f(x(i))\sum_i \bar w^{(i)} f(x^{(i)})는 0.9989 f(x(3))0.9989\,f(x^{(3)})에 weight 합이 0.0011인 나머지 네 particle의 항이 더해진 값입니다.

Weight degeneracy는 normalized weight가 소수의 particle에 몰려서, particle 수 NN에 비해 estimate에 실제로 기여하는 particle이 몇 개뿐인 상태입니다. 기여하는 particle이 적으면 estimate가 그 몇 개의 위치에 크게 좌우되므로 variance가 커집니다.

곱셈으로 쌓이는 차이

§1의 곱 ∏kp(yk∣xk(i))\prod_k p(y_k\mid x_k^{(i)})를 두 trajectory로 비교합니다. Trajectory AA는 매 step likelihood 0.8을, BB는 0.6을 받는다고 합시다. 한 step에서는 비가 0.8/0.6≈1.330.8/0.6\approx1.33으로 크지 않습니다. 10 step 뒤의 비는 다음과 같습니다.

w10Aw10B=(0.80.6)10≈17.8\frac{w_{10}^{A}}{w_{10}^{B}}=\left(\frac{0.8}{0.6}\right)^{10}\approx 17.8

두 trajectory만 있다면 normalized weight는 AA가 0.947, BB가 0.053입니다. 50 step 뒤에는 비가 약 1.8×1061.8\times10^6이 됩니다.

한 번 작아진 weight는 잘 회복되지 않습니다. 현재 wA=0.999999w^A=0.999999, wB=0.000001w^B=0.000001이고, 다음 measurement에서 BB가 AA보다 100배 잘 맞아 likelihood가 LA=0.01L^A=0.01, LB=1L^B=1이라고 합시다.

wA∝0.999999×0.01≈0.01,wB∝0.000001×1=0.000001w^A \propto 0.999999\times0.01\approx0.01,\qquad w^B\propto 0.000001\times1=0.000001

Normalize하면 AA가 0.9999, BB가 0.0001입니다. SIS의 weight는 trajectory 전체의 기록이므로, 이번 measurement 하나를 훨씬 잘 설명해도 과거에 쌓인 작은 weight를 뒤집지 못합니다.

몰리는 것 자체는 올바른 계산이다

여기서 한 가지를 구분해야 합니다. 매 step 데이터를 더 잘 설명한 trajectory가 path posterior p(x0:t∣y1:t)p(x_{0:t}\mid y_{1:t})에서 더 큰 확률을 갖는 것은 Bayes theorem이 말하는 그대로입니다. 위의 17.8이라는 비율은 계산 오류가 아닙니다.

문제는 그 posterior를 NN개의 proposal path로 표현한다는 데 있습니다. Path space의 차원은 state 차원 nn에 대해 n(t+1)n(t+1)이고 시간에 따라 커집니다. Path posterior는 measurement가 쌓일수록 좁은 영역에 모이는데, proposal path는 measurement를 보지 않고 transition model로만 뽑힙니다.

단순화한 가정으로 대략적인 크기를 계산해 봅니다. 매 step transition sample이 그 step의 measurement와 맞는 영역에 들어갈 확률이 0.5이고 step끼리 independent하다고 가정합니다. 그러면 path 1000개 중 tt step 모두 맞는 path 수의 expectation은 1000×0.5t1000\times0.5^t이고, t=10t=10이면 0.98개, t=20t=20이면 0.00095개입니다. t=20t=20에서는 posterior가 모인 영역에 path가 거의 하나도 없고, weight는 "가장 덜 틀린" 몇 개에 몰립니다.

정확히 말하면 SIS의 weight는 차원이 tt와 함께 커지는 space 위의 importance ratio p(x0:t∣y1:t)/qt(x0:t)p(x_{0:t}\mid y_{1:t})/q_t(x_{0:t})이고, proposal이 target과 같지 않으면 이 ratio의 variance가 일반적으로 시간에 따라 커집니다. 이 강의에서 증명하지는 않습니다. 0.8과 0.6의 곱이 벌어지는 현상은 이 variance 증가가 weight에 나타나는 모습입니다. Importance Sampling in Off-Policy RL에서 action sequence 위의 ratio가 step별 ratio의 곱이라 horizon이 길수록 variance가 커진 것과 같은 구조입니다.

Degeneracy를 한 문장으로

Weight가 데이터를 잘 설명한 trajectory에 몰리는 것은 path posterior를 그대로 따른 결과입니다. 문제는 차원이 계속 커지는 path space에서 고정된 NN개의 proposal path가 posterior가 몰린 곳에 거의 떨어지지 않는다는 점이고, 그래서 weight의 variance가 커지고, estimate에 실제로 기여하는 particle이 몇 개로 줄어듭니다.

한 step 안에서도 생긴다

시간에 따른 누적이 없어도 degeneracy가 생기는 경우가 있습니다. Predicted particle 1000개가 [0,100][0,100]에 고르게 퍼져 있고, sensor가 정확해서 yt=50y_t=50이 들어오고 likelihood가 N(yt; xt, 0.52)\mathcal N(y_t;\,x_t,\,0.5^2)라고 합시다. ∣xt−50∣=2|x_t-50|=2인 곳의 likelihood는 최댓값(xt=50x_t=50)의 e−8≈0.00034e^{-8}\approx0.00034배이므로, 50에서 2 이상 떨어진 particle의 weight는 xt=50x_t=50에 있는 particle의 0.00034배 이하입니다.

∣xt−50∣≤1|x_t-50|\le1 안에 있는 particle 수의 expectation은 1000×2/100=201000\times2/100=20개입니다. Seed를 고정해 한 번 뽑아 보니 16개였고, §3의 effective sample size는 15.9였습니다. 첫 step에서도, resampling 직후에도 이런 일은 똑같이 일어납니다.

원인은 likelihood가 transition이 퍼뜨린 폭보다 훨씬 날카롭다는 것, 즉 한 step 안의 proposal mismatch입니다. 시간에 따른 누적과는 다른 원인이고, §5에서 보듯 resampling이 고치는 것은 누적 쪽뿐입니다.


3. Effective sample size

정의

Degeneracy가 얼마나 심한지 숫자 하나로 보는 흔한 지표가 있습니다.

  Neff=1∑i=1N(wˉ(i))2  \boxed{\;N_{\text{eff}}=\frac{1}{\sum_{i=1}^N\big(\bar w^{(i)}\big)^2}\;}

"Normalized weight 제곱합의 역수"라고 읽습니다. 극단적인 경우부터 계산합니다.

  • 모든 weight가 같으면 wˉ(i)=1/N\bar w^{(i)}=1/N이고 ∑i(1/N)2=N⋅1/N2=1/N\sum_i(1/N)^2=N\cdot 1/N^2=1/N이므로 Neff=NN_{\text{eff}}=N입니다.
  • 한 particle만 weight 1이면 제곱합이 1이므로 Neff=1N_{\text{eff}}=1입니다.
  • kk개가 1/k1/k씩 나눠 갖고 나머지가 0이면 제곱합이 k⋅1/k2=1/kk\cdot1/k^2=1/k이므로 Neff=kN_{\text{eff}}=k입니다.

그래서 NeffN_{\text{eff}}는 "weight가 고른 particle 몇 개와 비슷한 상태인가"로 읽습니다. 값은 항상 1과 NN 사이에 있습니다.

∑i(wˉ(i))2≤∑iwˉ(i)=1(0≤wˉ(i)≤1)1=(∑iwˉ(i))2≤N∑i(wˉ(i))2(Cauchy–Schwarz)\begin{aligned} \sum_i \big(\bar w^{(i)}\big)^2 &\le \sum_i \bar w^{(i)} = 1 && \big(0\le\bar w^{(i)}\le1\big)\\ 1=\Big(\sum_i \bar w^{(i)}\Big)^2 &\le N\sum_i \big(\bar w^{(i)}\big)^2 && \text{(Cauchy–Schwarz)} \end{aligned}

첫 줄에서 Neff≥1N_{\text{eff}}\ge1, 둘째 줄에서 Neff≤NN_{\text{eff}}\le N이 나옵니다.

숫자로 보기

Normalized weight ∑i(wˉ(i))2\sum_i(\bar w^{(i)})^2 NeffN_{\text{eff}}
[0.2, 0.2, 0.2, 0.2, 0.2][0.2,\,0.2,\,0.2,\,0.2,\,0.2] 0.2 5
[0.05, 0.10, 0.60, 0.20, 0.05][0.05,\,0.10,\,0.60,\,0.20,\,0.05] 0.415 2.41
[0.01, 0.01, 0.95, 0.02, 0.01][0.01,\,0.01,\,0.95,\,0.02,\,0.01] 0.9032 1.11
[0.0001, 0.0002, 0.9989, 0.0005, 0.0003][0.0001,\,0.0002,\,0.9989,\,0.0005,\,0.0003] 0.9978 1.002

§2의 두 trajectory는 10 step 뒤 [0.947, 0.053][0.947,\,0.053]이므로 2개 중 Neff=1.11N_{\text{eff}}=1.11이고, [0,100][0,100] 예시는 1000개 중 15.9입니다.

NeffN_{\text{eff}}는 weight만으로 계산하는 근사 지표입니다. Target에서 직접 뽑은 i.i.d. sample 몇 개와 같은 정확도인지를 대략 나타내도록 만든 식이고, 그 수를 정확히 주지는 않습니다. Particle의 위치는 보지 않으므로, 모든 particle이 posterior가 큰 영역에서 벗어나 있어도 weight가 고르면 큰 값이 나올 수 있습니다.

Adaptive resampling

NeffN_{\text{eff}}가 threshold NthrN_{\text{thr}}보다 작을 때만 resampling하는 방법이 흔합니다. 흔히 Nthr=N/2N_{\text{thr}}=N/2를 씁니다. 한 번도 resampling하지 않으면 SIS이고, 매 step 하면 §7의 bootstrap filter입니다. Threshold는 두 문제 사이의 선택입니다. 너무 자주 resampling하면 §6의 diversity를 잃고, 너무 드물게 하면 degeneracy가 쌓여 mass가 거의 없는 particle에 계산을 씁니다.

Resampling 직후에는 모든 weight가 1/N1/N이므로 정의에 따라 Neff=NN_{\text{eff}}=N입니다. 이 값은 서로 다른 값이 몇 개인지 알려주지 않습니다. 5개가 모두 x(3)x^{(3)}의 복사본이어도 Neff=5N_{\text{eff}}=5이지만 서로 다른 값은 하나입니다. Resampling 뒤에는 서로 다른 particle의 개수를 따로 봐야 하고, §5의 그림 오른쪽이 그 값입니다.


4. Resampling

아이디어

Weighted particle 묶음 {x(i), wˉ(i)}i=1N\{x^{(i)},\,\bar w^{(i)}\}_{i=1}^N에서, 각 particle을 확률 wˉ(i)\bar w^{(i)}로 고르는 추첨을 NN번 해서 새 particle NN개를 만들고 weight를 모두 1/N1/N로 둡니다. 이것이 resampling입니다.

추첨은 index random variable로 씁니다.

I∼Categorical(wˉ(1),…,wˉ(N)),P(I=i)=wˉ(i),x~=x(I)I \sim \mathrm{Categorical}\big(\bar w^{(1)},\dots,\bar w^{(N)}\big),\qquad P(I=i)=\bar w^{(i)},\qquad \tilde x = x^{(I)}

x~\tilde x는 "resampling으로 만든 새 particle"이라고 읽습니다. 새 particle은 새로운 값이 아니라 기존 particle 중 하나의 복사본입니다.

Categorical sample 뽑기

Normalized weight가 wˉ=[0.05, 0.10, 0.60, 0.20, 0.05]\bar w=[0.05,\,0.10,\,0.60,\,0.20,\,0.05]라고 합시다. 누적합 Ci=∑j≤iwˉ(j)C_i=\sum_{j\le i}\bar w^{(j)}를 만들고 C0=0C_0=0으로 둡니다.

C=[0.05, 0.15, 0.75, 0.95, 1.00]C=[0.05,\,0.15,\,0.75,\,0.95,\,1.00]

U∼Uniform(0,1)U\sim\mathrm{Uniform}(0,1)을 하나 뽑고, UU 이상인 첫 누적합의 번호를 고릅니다.

I=min⁡{i:U≤Ci}I=\min\{i : U\le C_i\}

이 규칙에서 I=iI=i는 Ci−1<U≤CiC_{i-1}<U\le C_i와 같습니다. 그래서 구간은 (0, 0.05]→1(0,\,0.05]\to1, (0.05, 0.15]→2(0.05,\,0.15]\to2, (0.15, 0.75]→3(0.15,\,0.75]\to3, (0.75, 0.95]→4(0.75,\,0.95]\to4, (0.95, 1]→5(0.95,\,1]\to5입니다. 이 규칙이 올바른 확률을 주는지 확인합니다.

P(I=i)=P(Ci−1<U≤Ci)(definition of I)=Ci−Ci−1(U∼Uniform(0,1))=wˉ(i)(definition of Ci)\begin{aligned} P(I=i) &= P(C_{i-1}<U\le C_i) && \text{(definition of } I)\\ &= C_i-C_{i-1} && (U\sim\mathrm{Uniform}(0,1))\\ &= \bar w^{(i)} && \text{(definition of } C_i) \end{aligned}

U=0.63U=0.63이면 0.15<0.63≤0.750.15<0.63\le0.75이므로 I=3I=3입니다. U=0.82U=0.82이면 I=4I=4, U=0.03U=0.03이면 I=1I=1입니다. Weight가 큰 particle일수록 [0,1][0,1]에서 차지하는 구간이 넓으므로 더 자주 뽑힙니다.

bpf l7 resampling cdf
누적합 CiC_i로 [0,1][0,1]을 weight 길이만큼 나눈 띠입니다. UU가 떨어진 구간의 번호가 II이고, 세 화살표는 본문의 U=0.03, 0.63, 0.82U=0.03,\ 0.63,\ 0.82입니다.

다섯 번 뽑기

N=5N=5개를 다시 만들기 위해 UU를 다섯 번 independent하게 뽑았더니 0.63, 0.82, 0.41, 0.19, 0.120.63,\ 0.82,\ 0.41,\ 0.19,\ 0.12가 나왔습니다. Index는 3, 4, 3, 3, 23,\ 4,\ 3,\ 3,\ 2이고, 새 particle 묶음은 [x(3), x(4), x(3), x(3), x(2)][x^{(3)},\,x^{(4)},\,x^{(3)},\,x^{(3)},\,x^{(2)}]이며 weight는 모두 1/5=0.21/5=0.2입니다.

ii wˉ(i)\bar w^{(i)} 복제 수 NiN_i 새 묶음에서의 mass Ni/NN_i/N
1 0.05 0 0
2 0.10 1 0.2
3 0.60 3 0.6
4 0.20 1 0.2
5 0.05 0 0

x(3)x^{(3)}은 3개로 복제되어 mass가 3×0.2=0.63\times0.2=0.6이고 원래 weight와 같습니다. 이것은 이 realization에서 우연히 맞은 것입니다. 같은 표에서 x(2)x^{(2)}는 weight 0.10인데 mass 0.2를 받았고, x(1)x^{(1)}과 x(5)x^{(5)}는 0.05씩 잃었습니다. x(3)x^{(3)}이 정확히 3개 나올 확률도 (53)0.63 0.42≈0.346\binom{5}{3}0.6^3\,0.4^2\approx0.346뿐입니다.

복제 수는 random이다

NN번의 추첨이 independent하므로, 복제 수 (N1,…,NN)(N_1,\dots,N_N)은 Multinomial(N; wˉ(1),…,wˉ(N))\mathrm{Multinomial}(N;\,\bar w^{(1)},\dots,\bar w^{(N)})을 따릅니다. 이 방식을 multinomial resampling이라고 합니다. Ni=∑k=1N1{Ik=i}N_i=\sum_{k=1}^N\mathbf 1\{I_k=i\}로 쓰고 expectation을 계산합니다.

E[Ni]=∑k=1NE[1{Ik=i}](linearity of expectation)=∑k=1NP(Ik=i)(expectation of an indicator)=N wˉ(i)(P(Ik=i)=wˉ(i))\begin{aligned} \mathbb E[N_i] &= \sum_{k=1}^N \mathbb E\big[\mathbf 1\{I_k=i\}\big] && \text{(linearity of expectation)}\\ &= \sum_{k=1}^N P(I_k=i) && \text{(expectation of an indicator)}\\ &= N\,\bar w^{(i)} && \big(P(I_k=i)=\bar w^{(i)}\big) \end{aligned}
  E[Ni]=N wˉ(i),Var(Ni)=N wˉ(i)(1−wˉ(i))  \boxed{\;\mathbb E[N_i]=N\,\bar w^{(i)},\qquad \mathrm{Var}(N_i)=N\,\bar w^{(i)}\big(1-\bar w^{(i)}\big)\;}

Variance는 NiN_i가 성공 확률 wˉ(i)\bar w^{(i)}인 시행 NN번의 성공 횟수, 즉 binomial이라서 나옵니다. 같은 다섯 점에서 1000개를 뽑으면 expectation은 50, 100, 600, 200, 50개이고 standard deviation은 6.9, 9.5, 15.5, 12.6, 6.9개입니다. Seed를 고정해 한 번 뽑아 보니 50, 96, 588, 208, 58개였습니다.

Resampling이 weight를 복제 수로 바꾼다는 말은 expectation에서 성립합니다. 한 번의 결과 NiN_i는 expectation과 다를 수 있고, 그 퍼짐의 크기가 위의 Var(Ni)\mathrm{Var}(N_i)입니다.

Resampling은 무엇에서 뽑는가

새 particle 하나가 x(i)x^{(i)}가 될 확률은 wˉ(i)\bar w^{(i)}입니다. 이것은 L6 §5의 weighted empirical distribution p^(x)=∑iwˉ(i)δ(x−x(i))\hat p(x)=\sum_i\bar w^{(i)}\delta(x-x^{(i)})에서 sample 하나를 뽑는 것과 같습니다.

Resampling은 진짜 posterior가 아니라 weighted particle 묶음이 만든 근사 distribution에서 뽑습니다. 그래서 resampled 묶음은 원래 묶음에 없던 값을 만들지 못합니다. 현재 시각의 expectation을 추정하는 데에도 원래 weighted 묶음보다 나아질 수 없으며, 이것은 §5에서 variance로 확인합니다. 대신 평균적으로는 원래 weighted estimate를 보존합니다.

E[1N∑k=1Nf(x~(k))  ∣  {x(i),wˉ(i)}]=1N∑k=1NE[f(x(Ik))  ∣  {x(i),wˉ(i)}](linearity of expectation)=1N∑k=1N∑i=1Nwˉ(i)f(x(i))(P(Ik=i)=wˉ(i))=∑i=1Nwˉ(i)f(x(i))\begin{aligned} \mathbb E\Big[\frac1N\sum_{k=1}^N f(\tilde x^{(k)})\;\Big|\;\{x^{(i)},\bar w^{(i)}\}\Big] &= \frac1N\sum_{k=1}^N \mathbb E\big[f(x^{(I_k)})\;\big|\;\{x^{(i)},\bar w^{(i)}\}\big] && \text{(linearity of expectation)}\\ &= \frac1N\sum_{k=1}^N \sum_{i=1}^N \bar w^{(i)} f(x^{(i)}) && \big(P(I_k=i)=\bar w^{(i)}\big)\\ &= \sum_{i=1}^N \bar w^{(i)} f(x^{(i)}) \end{aligned}

조건에 weighted 묶음이 들어가 있다는 점이 중요합니다. Weighted 묶음을 고정했을 때 unbiased라는 뜻이고, 그 묶음 자체의 오차는 그대로 남습니다.

Systematic resampling

실제 구현에서는 UU를 NN번 따로 뽑지 않는 방법도 많이 씁니다. Systematic resampling은 U1∼Uniform(0,1/N)U_1\sim\mathrm{Uniform}(0,1/N) 하나만 뽑고 Uk=U1+(k−1)/NU_k=U_1+(k-1)/N으로 간격이 1/N1/N인 점 NN개를 만든 뒤, 같은 규칙 Ik=min⁡{i:Uk≤Ci}I_k=\min\{i:U_k\le C_i\}을 씁니다. 길이 wˉ(i)\bar w^{(i)}인 구간에 간격 1/N1/N인 점이 들어가는 개수는 Nwˉ(i)N\bar w^{(i)}의 내림이나 올림뿐이므로, E[Ni]=Nwˉ(i)\mathbb E[N_i]=N\bar w^{(i)}은 그대로이면서 복제 수 NiN_i의 variance가 multinomial보다 크지 않습니다.

위 예시에서 U1=0.07U_1=0.07이면 점은 0.07, 0.27, 0.47, 0.67, 0.870.07,\ 0.27,\ 0.47,\ 0.67,\ 0.87이고 index는 2, 3, 3, 3, 42,\ 3,\ 3,\ 3,\ 4입니다. Nwˉ(3)=3N\bar w^{(3)}=3, Nwˉ(4)=1N\bar w^{(4)}=1이므로 x(3)x^{(3)}은 어떤 U1U_1에서도 정확히 3개, x(4)x^{(4)}는 정확히 1개입니다.


5. Resampling이 고치는 것과 대가

고치는 것: 시간에 쌓인 weight

Resampling 뒤에는 모든 particle의 weight가 1/N1/N이므로, 다음 step의 weight는 §1의 식에서 wt−1(i)w_{t-1}^{(i)}이 빠진 p(yt∣xt(i))p(y_t\mid x_t^{(i)})뿐입니다. §2에서 본 "과거의 작은 weight를 계속 곱해 들고 가는" 구조가 끊깁니다. Weight가 거의 0이던 particle은 대부분 복제 수 0으로 사라지고, weight가 큰 particle은 여러 복사본이 되어 각자 다음 transition을 뽑습니다.

이후 weight는 현재 particle들의 likelihood만으로 정해지고, 과거 weight의 차이는 weight에 남지 않고 복제 수에 반영됩니다. 그래서 SIS에서 시간에 따라 줄어들던 weighting 직후의 NeffN_{\text{eff}}가, 매 step resampling하면 step이 지나도 누적해서 줄어들지 않습니다.

같은 N, 같은 비용, 다른 배치

Weight가 [0.01, 0.01, 0.95, 0.02, 0.01][0.01,\,0.01,\,0.95,\,0.02,\,0.01]인 particle 5개를 생각합니다. Resampling하지 않으면 다음 step에서도 5개를 모두 propagate하고 likelihood를 계산하는데, 그중 4개는 합쳐서 mass 0.05만 가진 가설입니다. Resampling하면 대부분 x(3)x^{(3)}의 복사본이 되고, 다음 transition에서 각자 다른 noise를 받아 x(3)x^{(3)} 주변의 서로 다른 점에 놓입니다.

Particle 수는 resampling 전후 모두 NN이고, propagate와 weight의 비용도 step마다 O(N)O(N) 그대로입니다. Resampling 자체도 계산이 들며, systematic resampling은 O(N)O(N)입니다. Resampling은 계산량을 줄이지 않고, 같은 NN개의 계산을 posterior mass가 있는 곳에 다시 배치합니다. 그 이득은 다음 step부터 나타납니다.

고치지 못하는 것: 한 step 안의 mismatch

§2의 [0,100][0,100] 예시를 매 step resampling하는 filter에 넣어도, 그 step에서 propagate된 particle에 날카로운 likelihood를 곱하는 순간 weight는 다시 소수에 몰립니다. Resampling은 weighting이 끝난 뒤에 일어나므로 이 step의 mismatch를 바꾸지 못합니다.

도움이 되는 것은 셋입니다.

  1. yty_t를 보고 particle을 뽑는 proposal입니다. 이 시리즈의 범위 밖입니다.
  2. Particle 수를 늘리는 것입니다. 다만 likelihood가 날카로울수록 필요한 NN이 커지고, likelihood가 degenerate하면 어떤 NN으로도 해결되지 않습니다. 예를 들어 yt=xty_t=x_t처럼 sensor noise가 0이고 transition이 continuous density이면, particle이 정확히 xt=ytx_t=y_t에 떨어질 확률이 0이므로 probability 1로 모든 weight가 0이 됩니다.
  3. Measurement noise를 실제보다 크게 잡은 likelihood를 쓰는 것입니다. Weight가 덜 몰려서 실제로 estimate가 좋아지는 경우가 많지만, §6의 roughening처럼 filter가 근사하는 posterior가 원래 model의 posterior에서 달라집니다.

아래 그림은 1차원 random walk에 N=500N=500 particle을 쓴 결과입니다. 왼쪽에서 t=1t=1에는 두 방법이 같은 알고리즘이므로 NeffN_{\text{eff}}도 둘 다 227이고, 그 뒤 SIS는 t=6t=6에 2 아래로 떨어집니다. 매 step multinomial resampling하는 SIR에서는 weighting 직후의 NeffN_{\text{eff}}가 N=500N=500으로 돌아가지 않고 평균 199 (범위 31–284)에 머무는데, 이것은 이 한 step 안의 mismatch 때문입니다.

대가: 지금 step의 Monte Carlo variance

§4에서 resampled estimate의 conditional expectation은 weighted estimate와 같았습니다. 그러나 추첨의 randomness가 variance를 더합니다. Multinomial resampling이면 x~(k)\tilde x^{(k)}들은 weighted 묶음이 주어졌을 때 p^\hat p에서 뽑은 i.i.d. sample이므로, L5 §4의 계산을 그대로 씁니다.

Var(1N∑k=1Nf(x~(k))  ∣  {x(i),wˉ(i)})=1N Varp^(f(x~))(multinomial resampling: conditionally i.i.d. draws)=1N(∑iwˉ(i)f(x(i))2−(∑iwˉ(i)f(x(i)))2)(variance under p^)\begin{aligned} \mathrm{Var}\Big(\frac1N\sum_{k=1}^N f(\tilde x^{(k)})\;\Big|\;\{x^{(i)},\bar w^{(i)}\}\Big) &= \frac1N\,\mathrm{Var}_{\hat p}\big(f(\tilde x)\big) && \text{(multinomial resampling: conditionally i.i.d. draws)}\\ &= \frac1N\Big(\sum_i \bar w^{(i)} f(x^{(i)})^2-\Big(\sum_i \bar w^{(i)} f(x^{(i)})\Big)^2\Big) && \text{(variance under } \hat p) \end{aligned}

Systematic resampling에서는 x~(k)\tilde x^{(k)}들이 conditionally i.i.d.가 아니므로 conditional variance가 이 식과 다르지만, 0 이상이라는 점은 같습니다.

Conditional expectation이 weighted estimate와 같은 resampling 방식이면, law of total variance에 따라 resampled estimate의 전체 variance는 weighted estimate의 variance에 conditional variance의 expectation을 더한 값입니다. 그러므로 resampled estimate의 variance는 weighted estimate보다 작아질 수 없습니다. 그래서 시각 tt의 state estimate는 resampling 전의 weighted 묶음으로 ∑iwˉt(i)f(xt(i))\sum_i\bar w_t^{(i)}f(x_t^{(i)})를 계산합니다. Bias and Variance의 용어로는, resampling은 bias를 더하지 않고 현재 step의 variance를 더합니다.

추첨의 randomness가 만드는 대가는 variance만이 아닙니다. 추첨이 운 나쁘게 나오면 참 state 근처의 particle이 한 번에 모두 사라질 수 있고, 이것을 particle deprivation이라고 합니다. 한 번의 resampling에서 이 확률은 보통 NN이 커질수록 빠르게 작아지지만 매번 0보다 크므로, 오래 돌릴수록 한 번이라도 일어날 확률이 커집니다. Particle 수가 belief가 퍼진 영역을 덮기에 부족할수록 잘 일어납니다.

Resampling 뒤에 random particle을 몇 개 섞어 넣으면 deprivation을 줄일 수 있지만 완전히 막지는 못하고, roughening처럼 filter가 근사하는 posterior를 바꿉니다. Impoverishment는 올바른 영역에 같은 값이 중복된 상태이고, deprivation은 올바른 영역에 particle이 아예 없는 상태입니다.

또 하나의 대가는 서로 다른 particle 값이 줄어드는 것이고, §6에서 다룹니다.

Resampling은 degeneracy를 없애는 연산이 아니다

Resampling이 없애는 것은 시간에 따라 곱해진 weight입니다. Likelihood가 transition의 폭보다 날카로워서 생기는 한 step 안의 degeneracy는 그대로 남고, 이것은 더 나은 proposal, 더 많은 particle, 실제보다 크게 잡은 measurement noise로 줄여야 하며, likelihood가 degenerate하면 particle 수로는 해결되지 않습니다. 또 resampling은 현재 step의 variance를 더하고 distinct particle을 줄이며, 직후의 Neff=NN_{\text{eff}}=N은 weight가 같다는 것만 나타내고 distinct particle 수와는 무관합니다.

bpf l7 neff
1차원 random walk xt=xt−1+wtx_t=x_{t-1}+w_t, wt∼N(0,1)w_t\sim\mathcal N(0,1), yt=xt+vty_t=x_t+v_t, vt∼N(0,0.52)v_t\sim\mathcal N(0,0.5^2), x0∼N(0,1)x_0\sim\mathcal N(0,1)에 N=500N=500 particle을 seed 고정으로 돌린 결과입니다. 왼쪽은 SIS와 SIR의 weighting 직후 NeffN_{\text{eff}}이고 두 run은 첫 resampling 전까지 같은 random number를 쓰며, 오른쪽은 SIR에서 resampling 뒤 서로 다른 particle 수입니다.

오른쪽에서 resampling 뒤 서로 다른 particle은 평균 약 197개입니다. 같은 시점의 NeffN_{\text{eff}}는 정의상 500이지만, 서로 다른 값은 그보다 훨씬 적습니다 (§3).


6. Sample impoverishment와 diversity

복제가 만드는 문제

Weight가 [0.01, 0.01, 0.96, 0.01, 0.01][0.01,\,0.01,\,0.96,\,0.01,\,0.01]이면 multinomial resampling의 결과가 [x(3), x(3), x(3), x(3), x(3)][x^{(3)},\,x^{(3)},\,x^{(3)},\,x^{(3)},\,x^{(3)}]일 확률이 0.965≈0.8150.96^5\approx0.815입니다. Particle은 5개이고 weight도 고르지만 값은 하나뿐입니다.

Sample impoverishment는 resampling이 반복되면서 particle이 가진 서로 다른 값의 수가 줄어드는 현상입니다. Weight degeneracy가 "값은 다양한데 weight가 몰린 상태"라면, impoverishment는 "weight는 고른데 값이 중복된 상태"입니다. Degeneracy가 심할수록 resampling이 소수를 많이 복제하므로 impoverishment도 심해집니다.

Transition noise가 다시 퍼뜨린다

Resampling 뒤 세 particle이 모두 xt=5.0x_t=5.0이고, transition model이 xt+1=xt+1+ϵx_{t+1}=x_t+1+\epsilon, ϵ∼N(0, 0.22)\epsilon\sim\mathcal N(0,\,0.2^2)라고 합시다. Particle마다 noise를 independent하게 뽑았더니 −0.10, 0.04, 0.18-0.10,\ 0.04,\ 0.18이 나왔습니다.

xt+1(1)=5.90,xt+1(2)=6.04,xt+1(3)=6.18x_{t+1}^{(1)}=5.90,\qquad x_{t+1}^{(2)}=6.04,\qquad x_{t+1}^{(3)}=6.18

같은 조상에서 출발했지만 다음 시각에는 세 개의 서로 다른 가설이 됩니다. 이 회복은 process noise가 degenerate하지 않을 때만 일어나고, 얼마나 퍼지는지는 noise의 크기가 posterior의 폭에 비해 얼마인지에 달려 있습니다. Noise가 posterior 폭보다 훨씬 작으면 복사본들은 거의 같은 자리에 머뭅니다.

퍼뜨릴 noise가 없을 때

모르는 sensor bias나 질량처럼 시간에 따라 변하지 않는 값을 state로 두면 transition model은 xt=xt−1x_t=x_{t-1}입니다. 이때 복사본은 다음 시각에도 똑같은 값이고, resampling을 할 때마다 서로 다른 값의 수는 줄어들 수만 있고 늘어나지 않습니다. Multinomial resampling을 반복하면 결국 probability 1로 모든 particle이 한 값의 복사본이 됩니다. 서로 다른 값이 kk개이면 가장 큰 mass를 가진 값의 mass가 1/k1/k 이상이므로, NN번의 추첨이 모두 그 값을 고를 확률은 k−N≥N−Nk^{-N}\ge N^{-N} 이상입니다. 값이 하나가 된 상태는 다시 벗어나지 못하고 매 resampling마다 이 양수 이상의 확률로 그 상태에 들어가므로, 반복하면 언젠가 들어갑니다.

Measurement가 없어 weight가 계속 같아도 이 일은 일어납니다. 그러면 filter는 한 값을 확신하는 것처럼 보이지만, 그 값은 아무것도 측정하지 않고 추첨으로 정해진 것입니다. Noise가 매우 작은 거의 deterministic한 dynamics도 비슷합니다.

먼저 resampling을 덜 할 수 있습니다. §3의 adaptive resampling은 weight가 고른 동안에는 resampling하지 않고 weight를 곱해서 들고 갑니다. §4의 systematic resampling은 weight가 모두 같으면 모든 particle을 정확히 한 번씩 남기므로, 이 경우 값을 잃지 않습니다.

그래도 resampling이 필요하면 resampling 뒤에 인공 noise를 더하는 방법이 있고, roughening이라고 합니다.

x(k)←x(k)+ϵrough(k),ϵrough(k)∼N(0, Σrough)x^{(k)} \leftarrow x^{(k)}+\epsilon_{\text{rough}}^{(k)},\qquad \epsilon_{\text{rough}}^{(k)}\sim\mathcal N(0,\,\Sigma_{\text{rough}})

Roughening은 model에 없는 noise를 transition에 더하는 것이므로, filter가 근사하는 posterior가 원래 model의 posterior에서 달라집니다. Target을 보존하면서 복사본을 흩뜨리는 MCMC move라는 방법도 있지만 이 시리즈에서는 다루지 않습니다.

Discrete state에서

State가 discrete이면 여러 particle이 같은 state에 있는 것 자체는 결함이 아닙니다. 그 state의 mass를 개수로 나타내는 방식입니다. Transition은 복사본을 한 step에 갈 수 있는 state로만, kernel의 확률대로 옮깁니다. 제자리에 머물 확률이 높으면 복사본 대부분이 같은 state에 남습니다.

Path diversity

시각 tt의 값이 다시 퍼져도, 각 particle의 조상을 거슬러 올라가면 resampling마다 조상 수가 줄어서 먼 과거에는 소수의 조상으로 모입니다. 이것은 trajectory x0:tx_{0:t} 전체를 추정하는 smoothing에서 문제가 되고, filtering distribution p(xt∣y1:t)p(x_t\mid y_{1:t})에는 직접 영향을 주지 않습니다.

두 문제 비교

Weight degeneracy Sample impoverishment
무엇인가 값은 다양한데 normalized weight가 소수에 몰림 weight는 고른데 값이 중복되어 서로 다른 값이 적음
원인 SIS에서 path space weight의 variance가 시간에 따라 커짐. 한 step 안에서는 likelihood가 transition 폭보다 날카로움 resampling의 복제, degeneracy가 심할수록 더 심함
도움이 되는 것 시간 누적에는 resampling. 한 step mismatch에는 yty_t를 쓰는 proposal, 더 많은 particle, measurement noise를 크게 잡은 likelihood resampling 빈도 줄이기 (adaptive resampling), systematic resampling, non-degenerate transition noise, roughening, MCMC move
부작용과 한계 resampling은 현재 step variance를 더하고 impoverishment와 deprivation 위험을 만듦. degenerate likelihood는 particle 수로 해결되지 않음. noise를 크게 잡으면 target이 바뀜 static parameter에는 transition이 도움이 안 됨. roughening은 target을 바꿈

7. SIR 알고리즘

한 cycle

지금까지의 단계를 한 시각 tt에 대해 순서대로 놓습니다.

  1. Propagate: xt(i)∼p(xt∣xt−1(i))x_t^{(i)}\sim p(x_t\mid x_{t-1}^{(i)}). Prediction 적분의 Monte Carlo 버전입니다 (L5 §6).
  2. Weight: wt(i)=wˉt−1(i) p(yt∣xt(i))w_t^{(i)}=\bar w_{t-1}^{(i)}\,p(y_t\mid x_t^{(i)}). 직전에 resampling했다면 wˉt−1(i)=1/N\bar w_{t-1}^{(i)}=1/N입니다 (§1).
  3. Normalize: wˉt(i)=wt(i)/∑jwt(j)\bar w_t^{(i)}=w_t^{(i)}/\sum_j w_t^{(j)}.
  4. Estimate: ∑iwˉt(i)f(xt(i))\sum_i\bar w_t^{(i)}f(x_t^{(i)}). Resampling 전에 계산합니다 (§5).
  5. Resample: 매 step, 또는 Neff<NthrN_{\text{eff}}<N_{\text{thr}}일 때 (§3, §4).
input : N
        sampler for p(x_0)
        sampler for p(x_t | x_{t-1})
        likelihood p(y_t | x_t)
        # N_thr > N: resample every step
        threshold N_thr
for i = 1..N:
    x_prev[i] ~ p(x_0)
    w_prev[i] = 1/N
for t = 1, 2, ...:
    receive y_t
    for i = 1..N:
        # 1. propagate
        x[i] ~ p(x_t | x_prev[i])
        # 2. weight
        w[i] = w_prev[i] * p(y_t | x[i])
    # s estimates p(y_t | y_{1:t-1})
    s = sum over i of w[i]
    for i = 1..N:
        # 3. normalize
        w_bar[i] = w[i] / s
    # 4. estimate, before resampling
    output sum over i of w_bar[i] * f(x[i])
    N_eff = 1 / (sum over i of w_bar[i]^2)
    # 5. resample
    if N_eff < N_thr:
        for k = 1..N:
            I ~ Categorical(w_bar[1..N])
            x_prev[k] = x[I]
            w_prev[k] = 1/N
    else:
        x_prev = x
        w_prev = w_bar

Neff≤NN_{\text{eff}}\le N이므로 Nthr>NN_{\text{thr}}>N으로 두면 매 step resampling합니다. 코드의 s는 버리는 값이 아닙니다. 매 step resampling하면 s는 1N∑ip(yt∣xt(i))\frac1N\sum_i p(y_t\mid x_t^{(i)})이고, 이 값은 L4 §5의 η−1=p(yt∣y1:t−1)\eta^{-1}=p(y_t\mid y_{1:t-1})을 추정합니다. L6 §4에서 self-normalization이 적분 ZZ를 피한 것이 아니라 같은 sample로 추정한 것과 같은 구조입니다.

이름

용어는 문헌마다 조금씩 다릅니다.

  • Transition model을 proposal로 쓰고 매 step resampling하는 filter를 bootstrap filter라고 합니다 (Gordon et al., 1993). Arulampalam et al. (2002)은 정확히 이 알고리즘을 SIR filter라고 부릅니다.
  • Doucet & Johansen (2009) 같은 문헌은 proposal이 무엇이든 resampling이 들어가면 SIR 또는 SISR이라고 부릅니다.
  • Resampling을 하지 않으면 SIS입니다.

약자 SIR은 sampling importance resampling으로도, sequential importance resampling으로도 풀어 씁니다. 이 시리즈에서 SIR은 따로 말하지 않는 한 bootstrap filter를 뜻합니다.

한 cycle을 그림으로

bpf l7 sir cycle
1차원 복도에서 SIR 한 cycle을 그린 것으로, predicted particle 20개를 두 후보 위치 3과 7 근처에서 seed 고정으로 뽑았습니다. Measurement model은 yt=xt+νty_t=x_t+\nu_t, νt∼N(0,0.52)\nu_t\sim\mathcal N(0,0.5^2)이고 받은 값은 yt=7.3y_t=7.3이며, 다음 transition은 xt+1=xt+1+ϵx_{t+1}=x_t+1+\epsilon, ϵ∼N(0,0.32)\epsilon\sim\mathcal N(0,0.3^2)입니다.

(2)에서 위치는 그대로이고 높이만 바뀌며, weighting 뒤 Neff=8.7N_{\text{eff}}=8.7입니다. (3)에서 왼쪽 무리는 복제 수가 모두 0이 되고 오른쪽 몇 개가 여러 번 복제되어, resampling 뒤 서로 다른 값은 20개 중 9개입니다. (4)에서 복사본이 independent한 noise를 받아 다시 서로 다른 값이 됩니다. 초록은 weight가 반영된 묶음, 파랑은 predicted particle이고, (4) 위 줄에서 같은 값의 복사본은 한 점으로 겹치므로 개수는 (3)의 높이로 읽습니다.

왼쪽 무리가 사라진 것은 measurement 7.3에서 그 가설들의 likelihood가 매우 작아 normalized weight 합이 약 2×10−142\times10^{-14}이기 때문이고, resampling은 그 결과를 particle 배치에 옮긴 것뿐입니다.


8. Bayes filter와 particle filter 대응

단계 Bayes filter Particle filter (SIR) 근사가 들어가는 곳
Belief 함수 belt−1(xt−1)\mathrm{bel}_{t-1}(x_{t-1}) 묶음 {xt−1(i), wˉt−1(i)}\{x_{t-1}^{(i)},\,\bar w_{t-1}^{(i)}\}, 즉 ∑iwˉt−1(i) δ(xt−1−xt−1(i))\sum_i \bar w_{t-1}^{(i)}\,\delta(x_{t-1}-x_{t-1}^{(i)}) 유한한 NN
Prediction bel‾t(xt)=∫p(xt∣xt−1) belt−1(xt−1) dxt−1\overline{\mathrm{bel}}_t(x_t)=\int p(x_t\mid x_{t-1})\,\mathrm{bel}_{t-1}(x_{t-1})\,dx_{t-1} xt(i)∼p(xt∣xt−1(i))x_t^{(i)}\sim p(x_t\mid x_{t-1}^{(i)}) 적분 대신 mixture에서 sampling
Update belt(xt)∝p(yt∣xt) bel‾t(xt)\mathrm{bel}_t(x_t)\propto p(y_t\mid x_t)\,\overline{\mathrm{bel}}_t(x_t) wt(i)=wˉt−1(i) p(yt∣xt(i))w_t^{(i)}= \bar w_{t-1}^{(i)}\,p(y_t\mid x_t^{(i)}) 위치는 그대로, weight만 곱함
Normalization η−1=∫p(yt∣xt) bel‾t(xt) dxt\eta^{-1}=\int p(y_t\mid x_t)\,\overline{\mathrm{bel}}_t(x_t)\,dx_t wˉt(i)=wt(i)/∑jwt(j)\bar w_t^{(i)}=w_t^{(i)}/\sum_j w_t^{(j)}, 그리고 ∑jwt(j)≈η−1\sum_j w_t^{(j)}\approx\eta^{-1} self-normalization
Estimate ∫f(xt) belt(xt) dxt\int f(x_t)\,\mathrm{bel}_t(x_t)\,dx_t ∑iwˉt(i)f(xt(i))\sum_i \bar w_t^{(i)} f(x_t^{(i)}) Monte Carlo error
대응 없음 (없음) Resampling exact filter에는 필요 없는 단계

앞의 다섯 줄은 L4의 식을 particle 위에서 한 줄씩 근사한 것입니다. 마지막 줄만 exact Bayes filter에 대응하는 단계가 없습니다. Exact filter는 belief 함수 전체를 들고 가므로 weight가 몰릴 대상이 없고, resampling은 belief를 유한한 개수의 weighted particle로 표현하기 때문에 생긴 단계입니다.

이 시리즈의 식이 논문에서 어떻게 쓰이는지는 Lee 2020 · Prerequisites — Lecture Notes › 5 · Curriculum: from Bayes filter to particle filter에서 이어집니다.


이해 확인

1. SIS의 식 wt(i)∝wt−1(i) p(yt∣xt(i))w_t^{(i)}\propto w_{t-1}^{(i)}\,p(y_t\mid x_t^{(i)})을 유도할 때 target을 filtering distribution이 아니라 path posterior로 잡는 이유를 설명하시오. 또 매 step resampling하면 wt−1(i)w_{t-1}^{(i)}을 지워도 되는 이유를 설명하시오.

2. Trajectory AA가 매 step likelihood 0.8, BB가 0.6을 받아 10 step 뒤 weight 비가 약 17.8이 되었다. 이 비율이 계산 오류가 아니라면, SIS의 degeneracy에서 정확히 무엇이 문제인지 설명하시오.

3. 매 step resampling하는 SIR에서도 weighting 직후의 NeffN_{\text{eff}}가 NN보다 훨씬 작을 수 있다. 그 이유와 도움이 되는 방법을 설명하시오. 또 resampling 직후의 Neff=NN_{\text{eff}}=N을 보고 안심하면 안 되는 이유를 설명하시오.

4. 시각 tt의 estimate를 resampling 뒤가 아니라 앞에서 계산하는 이유를 설명하시오. 또 시간에 따라 변하지 않는 parameter를 state로 두고 resampling을 반복하면 어떤 일이 생기는지 설명하시오.


답

1. Particle이 자기 조상만 따라가면 뽑힌 것은 trajectory이고 그 proposal density가 transition의 곱으로 쌓이므로, importance ratio를 trajectory 위에서 세워야 합니다. 또 path posterior는 model factor의 곱에 비례해 점마다 계산되지만, exact filtering distribution은 bel‾t\overline{\mathrm{bel}}_t라는 적분이 필요합니다. Particle mixture로 근사해도, weight가 서로 다른 채로 조상만 따라가면 proposal이 그 mixture가 아니라서 L6 §6의 약분이 성립하지 않습니다. 매 step resampling하면 resampled path들은 resampling 전 묶음이 만든 empirical path distribution에서 뽑은 sample이므로, 이 근사를 target으로 쓰면 weight가 모두 1/N1/N이고 공통 상수로 ∝\propto에 흡수됩니다. → §1

2. 데이터를 더 잘 설명한 trajectory가 path posterior에서 커지는 것은 올바른 결과입니다. 문제는 path space의 차원이 tt와 함께 커져서, measurement를 보지 않고 뽑은 고정된 NN개의 path가 posterior가 몰린 좁은 영역에 거의 떨어지지 않는다는 것입니다. 그래서 weight의 variance가 커지고 NeffN_{\text{eff}}가 시간에 따라 1 근처까지 줄어듭니다. → §2

3. Likelihood가 transition이 퍼뜨린 폭보다 날카로우면 그 step의 weighting에서 바로 weight가 몰립니다. Resampling은 weighting 뒤에 일어나므로 이것을 바꾸지 못합니다. yty_t를 쓰는 proposal, 더 많은 particle, 실제보다 크게 잡은 measurement noise가 도움이 되지만, likelihood가 degenerate하면 particle 수로는 해결되지 않습니다. Resampling 직후의 Neff=NN_{\text{eff}}=N은 weight가 모두 1/N1/N이라 정의상 나오는 값이고, 복사본이 몇 종류인지는 알려주지 않으므로 distinct particle 수를 봐야 합니다. → §5, §3

4. Resampled estimate는 weighted estimate와 conditional expectation이 같지만 추첨의 variance가 더해지므로 더 나을 수 없습니다. Static parameter는 transition이 xt=xt−1x_t=x_{t-1}이라 복사본을 퍼뜨릴 noise가 없고, resampling마다 서로 다른 값이 줄어들기만 하고, multinomial resampling을 반복하면 결국 probability 1로 한 값만 남습니다. Roughening으로 흩뜨릴 수 있지만 target이 바뀝니다. → §5, §6

lectureprobabilitybayes-filterparticle-filter