論文読みメモ： Nonparametric density estimation in compound Poisson process using convolution power estimators（その1）

以下の論文を読みます。

F. Comte, C. Duval, V. Genon-Catalot. Nonparametric density estimation in compound Poisson process using convolution power estimators. Metrika, Springer Verlag, 77 163–183, 2014.
https://hal.archives-ouvertes.fr/hal-00780300/document

※ キャラクターは架空のものです。解釈の誤りは筆者に帰属します。お気付きの点がありましたらご指摘ください。

次回：まだ

f:id:cookie-box:20190101155733p:plain:w60

$(\xi_j, j \geqq 0)$ を独立に同一の確率密度 $f$ にしたがう実数値確率変数とし、 $(N_t)$ を強度 $c>0$ のポアソン過程として、複合ポアソン過程 $(X_t, t \geqq 0)$ は以下のように表せるということです。

$\displaystyle X_t = \sum_{i=1}^{N_t} \xi_i \tag{1}$ これはつまり、時刻 $t$ までに隕石が庭に降ってくる回数 $N_t$ が強度 $c$ のポアソン過程にしたがい、個々の隕石の質量 $\xi_i$ は独立に確率密度 $f$ にしたがうというときに、時刻 $t$ までに庭に降ってくる隕石の質量の総和ですね。

f:id:cookie-box:20190101160814p:plain:w60

喩えが大災害だな…。

しかし、隕石がいつ庭に落ちてくるかずっと監視しているほど暇ではありませんので、一定時間間隔 $\Delta$ ごとに、その時点までに降ってきた隕石の総質量 $(X_{j \Delta}, j \geqq 0)$ が観測されるのみとしましょう。このときに隕石の質量がしたがう確率密度 $f$ を推測するのが今回のやりたいことですね。無論 $c$ も $f$ も未知です。

いや庭に隕石が振ってきたら気付くよね普通…。

複合ポアソン過程はレヴィ過程の特別な場合で、レヴィ密度が $cf(\cdot)$ の場合に相当します。なので、もし $c$ が既知ならば $f$ の推定はジャンプのみからなるレヴィ過程のレヴィ密度の推定に等しいと。これについてはこれまでに多くの論文で取り扱われているとあります。しかし、複合ポアソン過程に特化した推定手法も研究されていて、Buchmann and Grübel (2003) は初めて特に $f$ を推定する手法を考案したのですかね？実際、1回の観測間隔の増分 $X_\Delta$ の確率密度は以下のようになるということです。

$\mathbb{P}_{X_\Delta}(dx) = e^{-c\Delta} \delta_0 (dx) + (1 - e^{-c\Delta}) g_\Delta (x) dx \tag{2}$ ただし $g_\Delta$ は $X_\Delta$ が $0$ でない下での条件付き分布であって、
$\displaystyle g_\Delta = \sum_{m \geqq 1} \frac{e^{-c\Delta}}{1 - e^{-c\Delta}} \frac{(c \Delta)^m}{m!} f^{\ast m} \tag{3}$ $f^{\ast m}$ は $f$ の $m$ 次の畳み込み（合成積）であるということですが…なんでこうなるんですか？

まず時刻 $\Delta$ までに隕石の降ってくる回数 $m$ がしたがう確率分布は $P(N_\Delta = m) = e^{-c\Delta} (c \Delta)^m / m!$ だから、時刻 $\Delta$ までに1個も隕石が振ってこない確率は $P(N_\Delta = 0) = e^{-c\Delta}$ で、1個以上隕石が振ってくる確率は $P(N_\Delta \geqq 1) = 1 - e^{-c\Delta}$ だね。さらに $m \geqq 0$ 回隕石が振ってくるとき、その質量の和 $\xi_1 + \xi_2 + \cdots + \xi_m$ がしたがう密度関数は、 $f$ を $m$ 回畳み込んだ関数になる。…あたりをつかえば、 $(2)$ と $(3)$ が成り立つことが容易にわかるよ。

合成積（畳み込み）の意味と応用３つ | 高校数学の美しい物語

だから増分 $X_\Delta$ をたくさん観測しても $X_\Delta = 0$ のデータについては $f$ を知る上で「何の足しにもならない」…その $\Delta$ 間に隕石が降ってこなかったら隕石の質量がしたがう分布について何も得るものがないのは当然だよね。だから van Es et al. (2007) はゼロでない観測だけ抜き出したらしい。以下のようにすれば増分がゼロでないインデックス $S_i$ だけ取ってこれるね。
$S_1 = {\rm inf} \{j \geqq 1, X_{j \Delta} - X_{(j-1)\Delta} \neq 0 \} , \quad S_i = {\rm inf} \{j > S_{i-1}, X_{j \Delta} - X_{(j-1)\Delta} \neq 0 \} \tag{4}$
$Z_i = X_{S_i \Delta} - X_{(S_i - 1) \Delta} \tag{5}$ こうして $(S_i, Z_i) \, (i = 1, \cdots, n)$ を集めれば、 $Z_1, \cdots, Z_n$ は $X_\Delta$ が $0$ でない下での条件付き分布 $g_\Delta$ にしたがう。解くべき問題は「確率密度 $g_\Delta$ からの $n$ 個のサンプルから $f$ を推定する」にブレークダウンされたわけだ。

van Es et al. (2007) は $c$ が既知で $\Delta = 1$ のケースについて、 $f$ の特性関数と $g_\Delta$ の特性関数の関係を利用して $f$ の推定量を構築したんですね。Duval (2012a) はまた別の推定量を考えたようです。変換 $f \to g_\Delta := P_\Delta f$ は明示的に逆変換できるので $f = P_\Delta^{-1} g_\Delta$ でよいと。 $c\Delta < \log2$ を仮定すれば $P_\Delta^{-1}$ は以下の級数で与えられるそうです。

$\displaystyle P_{\Delta}^{-1} (g) = \sum_{m \geqq 1} \frac{(-1)^{m+1}}{m} \frac{(e^{c\Delta} - 1)^m}{c \Delta} g^{\ast m} \tag{6}$ 実際には $K+1$ 項目まででトランケートして以下を推定量とすると。
$\displaystyle f \cong \sum_{m \geqq 1}^{K+1} \frac{(-1)^{m+1}}{m} \frac{(e^{c\Delta} - 1)^m}{c \Delta} g_{\Delta}^{\ast m} \tag{7}$ $\Delta$ が小さい場合はこの近似は有効であるということですが…なぜこうなるんでしょう？

…合成積ってフーリエ変換すると、合成積する前の関数のフーリエ変換の積になるよね。そう考えると、 $(3)$ って両辺をフーリエ変換すれば右辺がエクスポネンシャルのマクローリン展開にみえる。そう考えて式変形すると以下までたどり着けるかな。

$\displaystyle \mathcal{F}[f] = \frac{1}{c\Delta} \log \left( 1 + \frac{1 - e^{-c\Delta}}{e^{-c\Delta}} \mathcal{F}[g_\Delta] \right)$ そしたらこれをもっかいマクローリン展開すればいい。 $\displaystyle \mathcal{F}[f] = \frac{1}{c\Delta} \sum_{m \geqq 1} \frac{(-1)^{m+1}}{m} (e^{c\Delta} - 1)^m \mathcal{F}[g_\Delta]^{m}$ これを逆フーリエ変換すれば $(6)$ を得るね。

マクローリン展開 log(1+x)

普通に Duval (2012a) の原論文もみられるから貼っておくね。27～28ページのところでやっぱりフーリエ変換しているよね。

[1203.3135] Adaptive wavelet estimation of a compound Poisson process

そういうわけで、 $g_{\Delta}^{\ast m}$ を $Z_1, \cdots, Z_n$ から推定すればいいことになったけど、これはこれでそう簡単でもないっぽいな。Duval は wavelet threshold estimator というのを構築したらしい。 $g_{\Delta}^{\ast m}$ は $g_{\Delta}$ にしたがう確率変数 $m$ 個の和がしたがう密度なんだから、 $m$ 個の観測の和をとって推定したのかな？これはこれで収束性はよいみたいだけど、 $Z_1, \cdots, Z_n$ から $m$ 個のサンプルの和をつくるのは計算コストがかかると。だからこの論文では $(7)$ を利用するけど $g_{\Delta}^{\ast m}$ の推定手法は変えるらしい。

なるほど…先に進むと、以降、 $g_\Delta, \,g_{\Delta}^{\ast m}$ を表記するとき $\Delta$ を省くとあります。そして、「確率密度 $g$ からの $n$ 個のサンプルから、 $g^{\ast m} \, (m \geqq 2)$ の『 $\sqrt{n}$ -consistent nonparametric estimator』がつくれることはよく知られている」とありますね。これはどういう意味ですか？私によく知られていないんですが。

特に最近 Chesneau et al. (2013) がシンプルな推定量を提案したとあって以下に続くのがその方法なんじゃないかな。 $g^\ast$ を $g$ のフーリエ変換、 $(g^\ast)^m$ を $g^{\ast m}$ のフーリエ変換とするよ（ややこしいな）。この合成フーリエ変換（？） $(g^\ast)^m$ に対し、経験合成フーリエ変換（？） $(\tilde{g}^\ast)^m$ を以下で定義する。

$\displaystyle \tilde{g}^\ast(t) = \frac{1}{n} \sum_{j=1}^n e^{it Z_j} \tag{9}$ 本当は $g^\ast = \int g(x) e^{itx} dx$ なのを経験分布上での平均にした版だね。そして $(9)$ の $m$ 乗を逆フーリエ変換して $g^{\ast m}$ の推定量を得る。このときカットオフを $\ell$ として積分区間を $[-\pi \ell, \pi \ell]$ とするみたいだね。 $\displaystyle \widehat{g_{\ell}^{\ast m}} (x) = \frac{1}{2\pi} \int_{-\pi \ell}^{\pi \ell} e^{-itx} \bigl( \tilde{g}^\ast (t) \bigr)^m dt \tag{10}$ これをそのまま $(7)$ に代入すればいい。 $\displaystyle \widetilde{f_{K, \ell}}(x) = \sum_{m \geqq 1}^{K+1} \frac{(-1)^{m+1}}{m} \frac{(e^{c\Delta} - 1)^m}{c \Delta} \widehat{g_{\ell}^{\ast m}} (x) \tag{11}$

おお、これで隕石の質量がしたがう分布 $f$ が手に入り…ませんね。副部長、いま隕石が到着する頻度の強度 $c$ は未知です。右辺に $c$ があっては駄目でしょう。

いやだからまだこれは $f$ の推定量 $\widehat{f_{K, \ell}}(x)$ じゃないよ。 $\displaystyle c_m(\Delta) = \frac{(e^{c\Delta} - 1)^m}{c \Delta}$ にも何か推定量 $\widehat{c_m(\Delta)}$ をつくるみたいだね。以降で、 $\widehat{f_{K, \ell}}$ のL2リスクと、データに応じた $\hat{\ell}$ の選択方法が示されるみたいだよ。2節で $\widehat{c_m(\Delta)}$ をつくってそのL2リスクを議論して、3節で Chesneau et al. (2013) の合成確率密度の推定を復習して、4節で $f$ の推定とL2リスクを議論して、5節でさまざまな分布に対してシミュレーションして、6節がまとめ、7節に定理の証明が集められているのかな。

え、ああ、まだイントロダクションだったんですね…普通に番号付きの数式がどんどん出てくるのでてっきりもう本論に入っているのかと…。それでようやく2節に進むと、この節は Proposition 2つからなっていますね。1つ目の Proposition の主張は、 $k \geqq 1$ に対して、以下が成り立つということですか？

$\mathbb{P}(S_1 = k, Z_1 \leqq x) = e^{-c(k-1) \Delta} (1 - e^{-c \Delta}) \mathbb{P} (X_\Delta \leqq x | X_\Delta \neq 0)$ この数式は日本語訳するならば、

　「初めてゼロでない質量を観測するインデックス $S_1$ が $k$ に等しく、かつそのとき観測された質量が $x$ 以下である確率」は、
　「k-1回目までの観測ですべて『隕石が1個も降ってこない」を引く確率」
　×「k回目の観測で『隕石が1個は降ってくる』を引く確率」
　×「その観測での増分がゼロでない条件の下で、その観測の増分 $Z_1$ が $x$ 以下である確率」

に等しいということですね。したがって、初めてゼロでない質量を観測するインデックス $S_1$ と、そのとき観測される質量 $Z_1$ は独立だと。そして $S_1$ はパラメータが $1 - e^{-c \Delta}$ の幾何分布にしたがうと。これはそうなりそうですね。2つ目の Proposition の主張は、 $c \in [c_0, c_1], \, c_0 > 0, \, c_1 \Delta \leqq \log(2)/2$ として、 $1$ 以上の $m$ に対して次の $(12)$ $(13)$ を定義すると、 $\displaystyle H_m(\xi) = \frac{1}{(\xi - 1)^m \log \frac{\xi}{\xi - 1}} \tag{12}$
$\displaystyle \widehat{c_m(\Delta)} = H_m(S_n / n) {\bf 1}_{ \left\{1 + \frac{1}{e^{2 c_1 \Delta} -1} \leqq \frac{S_n}{n} \leqq 1 + \frac{1}{e^{c_0 / (2\Delta)} -1} \right\} } \tag{13}$ 以下が成り立つということですね。 $C_m$ は $c_0, c_1, m$ の関数であるそうです。 $\displaystyle \mathbb{E} \left( \widehat{c_m(\Delta)} - c_m(\Delta) \right)^2 \leqq C_m \frac{\Delta^{2(m-1)}}{n} \tag{14}$ つまり $\widehat{c_m(\Delta)}$ は漸近的ではなく普通に $c_m(\Delta)$ との2乗誤差が上式で抑えられる推定量になっているということですが…この関数をいきなりみてもどうしてこうなったのやらですね…。

証明は 9 ページにあるね。まず $p(\Delta) = 1 - e^{-c\Delta}$ とおいている。これは $\Delta$ 間に「『隕石が1個は降ってくる』を引く確率」だね。これは、 $\displaystyle c \Delta =- \log \bigl( 1 - p(\Delta) \bigr) = \log \left( \frac{p(\Delta)^{-1}}{p(\Delta)^{-1} - 1} \right)$ の関係がある。そうすると、 $\displaystyle c_m(\Delta) = \frac{(e^{c\Delta}-1)^m}{c \Delta} = \frac{(e^{c\Delta}-1)^m}{\log \left( \frac{p(\Delta)^{-1}}{p(\Delta)^{-1} - 1} \right)} = \frac{1}{\bigl(p(\Delta)^{-1}-1 \bigr)^m \log \left( \frac{p(\Delta)^{-1}}{p(\Delta)^{-1} - 1} \right)} = H_m \bigl( p(\Delta)^{-1} \bigr)$ となる。 $h_m$ は $\Delta$ 間に「『隕石が1個は降ってくる』を引く確率」（の逆数）を求める量に変換してくれる関数だったってだけだね。そして、 $\Delta$ 間に「『隕石が1個は降ってくる』を引く確率」の最尤推定量は $n / S_n$ になる。 $S_n$ は「隕石が1個は降ってくる」を $n$ 回観測するのに費やした観測の回数だからね。だから $\widehat{c_m(\Delta)} = H_m(S_n / n)$ としたいけど、 $S_n / n$ は $1$ になりうるのでこれはできない、か。確かに観測する度に隕石が1個以上落ちているケースだと、隕石が落ちたインデックスから強度を逆算しようとしても「常に隕石が落ちてる」にしかならないね。 $S_n / n$ が $1$ より大きくないと困る。だから、真の強度が入っている区間を $[c_0, c_1]$ だろうということにしておいて、観測された $S_n / n$ が想定される範囲をずれていたら $\widehat{c_m(\Delta)}$ はゼロにしてしまうという回避策を入れたみたいだ。実際、 $c_m(\Delta)$ は強度が大きいとゼロに近づくし…ただ、 $\widehat{c_m(\Delta)} = 0$ だと $f$ の推定はできなさそうだな…。