2021.02.24
EM法は確かに計算が簡単なので、有用である。しかし、原理的にはまだ判らないので、文献を調べた。
最初にエイズのデータに使用した Becker 等の論文が入手できないので、そもそも EM法を数学的にまとめあげた Dempster et al. "Maximum Liklihood from Incomplete Data via the EM Algorithm", Journal of the Royal Statistical Society: Series B (Methodological), 39(1) 1-2(1977) https://doi.org/10.1111/j.2517-6161.1977.tb01600.x の Introduction だけ読んだ。x から y への（多→１）写像が定義されている。それぞれがある確率分布に従うが、パラメータは不明である。y が実際に観測されたとして、その事実から、分布パラメータ Φ を決めるという問題である。その原理は尤度の最大化である。簡単な場合は直接計算できるが、一般的には解（Φ）が不定となる。これを反復計算で行う。まず試行解を設定して、その試行解（パラメータ）から x を推定する（Estimation）。その推定された x を観測データと見做して、x,y から尤度を与えて、それを最大化（Maximization）するように、改良されたパラメータ Φ を決める。尤度の表現の中にあるパラメータの一部を旧パラメータとして固定し、残りのパラメータについて尤度を最大化する、という事を繰り返せば、尤度は確実に大きくなっていく。
これだけではどういう風に応用すればよいのかが良く判らないので、2003年の SARS に応用した Yip et al. "Reconstrucition of the Infection Curve for SARS Epidemic in Beijing, China Using a Back-Projection Method", Communications in Statistics - Simulation and Computation, 37(2), 425-433(2008), https://doi.org/10.1080/03610910701792562 の 第２節（方法論）だけを読んだ。観測されるデータは日毎の発症者数 Y(t) である。t は最大で T までとする。発症者数は期待値 μ(t) というパラメータを持つポワッソン分布に従うとする。他方、日毎の感染者数は期待値 λ(t) というパラメータを持つポワッソン分布に従うとする。感染してから d 日後の発症確率を f(d) とする。μ(t)＝∑d λ(t-d)f(d) である。
ここまでは判っていたが、問題は Dempster 論文の x である。これが感染者数だと思っていたので理解できなかったのである。実は、これは、t 日に感染して、そこから d 日後の発症した数 N(t,d) （確率変数）となっている。上記 μ(t) の和の中の各項である。そして、この N(t,d) という疑似観測値に対してポワッソン分布を考える。その期待値パラメータは λ(t)f(d) となる。ポワッソン分布は再生性があるから、Yt に対するポワッソン分布の期待値パラメータが μ(t) ＝∑d N(t-d,d)となる。この確率変数 N(t,d) のポワッソン分布について、疑似観測値 N(t,d) を得た時の期待値パラメータ λ(t) についての対数尤度は、ポワッソン分布の関数形から、
（式ａ） logL＝∑(t=1:T)∑(d= 0:T-t) {N(t,d)log(λ(t)f(d)) - λ(t)f(d)}
となる。確率変数が Y(t) ではなく、N(t,d) になっている。ただ、一番外側の t の範囲としては、t=-L+1:T とすべきだろう（Lは最長の潜伏期間）。そうすると、その内側の d の範囲としては、t+d≧1 という条件が加わる。
１．Estimation Step
さて、最初の Estimation Step では N(t,d) の「期待値」を推定する。時刻 t+d での発症者数観測値 Y(t+d) を与える期待値は ∑(i=1:t+d) λ(i)f(t+d-i) であるから、この和の中で、t に感染して t+d に発症する場合の比率を考えると、
比率計算から、
（式ｂ） Expected N(t,d)＝Y(t+d) [λ(t)f(d)/{∑(i=1:t+d) λ(i)f(t+d-i)}]
となる。
２．Maximization Step
次は上記の対数尤度を最大化するが、ここで、N(t,d) を観測値と見做したのだから、その表式の中にある λ(t) は旧値 λ(t)old に固定する。最大化するために、残りの λ(t)new で微分して 0 と置く。各 t 独立に 0 としなくてはならないから、
∑(d=0:T-t) {N(t,d)/λ(t)new - f(d)}＝0
となり、当然ながら、和の中の各項を 0 にすると何も変化しない。しかし、d についての和があることで λ(t) が変化する。N(t,d) を元に戻すと、
∑(d= 0:T-t) [{Y(t+d)f(d)/∑(i= 1:t+d)λ(i)old f(t+d-i)} λ(t)old/λ(t)new - f(d)]＝0
これから、最終的に、
（式ｃ） λ(t)new/λ(t)old＝∑(d=0:T-t) [{Y(t+d)f(d)/∑(i=1:t+d) λ(i)old f(t+d-i)}] / {∑(d=0:T-t) f(d)}
が得られる。
これは簡便で、なかなか使えそうである。下記論文のように、X線CT等にも使われているようである。昔やった インピーダンスCT にも使えるかもしれない。そういえば、調査資料の中に backprojection というのがあったような気もするが、結果はあまり冴えなかったと記憶している。。。
Silverman et al., "A Smoothed Em Approach to Indirect Estimation Problems, with Particular Reference to Stereology and Emission Tomography", Journal of the Royal Statistical Society: Series B (Methodological) 52(2) 271-303(1990), https://doi.org/10.1111/j.2517-6161.1990.tb01788.x
ところで、Y(t+d) と {∑(i=1:t+d) λ(i)f(t+d-i)}] は一致しないし、この方法はそもそも一致することを目指していない。その比（誤差）を f(d) という重み付きで平均したものを 0（比で言うと 1）にするような方法である。つまり、方法の中に、時間についての f(d) というフィルター関数が含まれている。だから、非線形最小自乗法で β を大きくした場合と同じような結果になる。f(d) の拡がりが β で仕切られる固有関数の最短波長の半分位になるのであろう。