« 読了：Rizzi et al. (2016) ヒストグラムから元の確率分布をノンパラメトリックに推定する方法コンテスト | メイン | 読了: Camerer, et al. (2018) 社会科学の有名な実験研究21本を追試してみたら、ああなんてこったい、結果は... »
2018年10月16日 (火)
仕事の都合で突然訪れた「ヒストグラムから確率分布をノンパラメトリック推定したい」祭り、時間がないのでそろそろ終了したいんだけど、Rizzi et al.(2016)が挙げていたpenalized composite linkモデルが良く理解できず、なんだか気持ち悪いので、大急ぎで目を通した。Rのungroupパッケージの元論文である。
Rizzi, S., Gampe, J., Eilers, P.H. (2015) Efficient Estimation of Smooth Distributions From Coarsely Grouped Data. American J. Epidemiology, 182(2), 138-147.
イントロ... [省略]
ヒストグラムの元データに出現しうる値の系列(たとえば年齢だったら, 0歳, 1歳, ...)を$a_1, \ldots, a_J$とし、その頻度の期待値を$\gamma_j$とする。なにを云ってるのかというと、値が$a_j$である確率を$p_j$とし、仮にサイズ$N$の標本があるならば、$\gamma_j = Np_j$である。
$a_j$の頻度は平均$\gamma_j$のポアソン分布に従う。我々が推定したいのは$\mathbf{\gamma} = (\gamma_1, \ldots, \gamma_J)^T$である。ここまではいいっすね。
問題は、我々が観察しているのが$a_j$それぞれの頻度ではなく、もっと粗いビンにおける頻度分布だという点である。
ビンを$I$個、観察された頻度を$Y_1, \ldots, Y_I$とする。これもポアソン分布に従う。各ビンでの平均を$\mu_i$として
$P(Y_i = y_i) = \mu_i^{y_i} \exp(-\mu_i) / y_i!$
である。
$\mu_i$は$\gamma_j$をいくつかまとめたものだ、と考えて、$\mathbf{\mu} = (\mu_1, \ldots, \mu_I)^T$について
$\mathbf{\mu} = C \mathbf{\gamma}$
としよう。$C$は$I \times J$の行列で、要素は1か0。ビン$I$がもとの$J$番目の値を含んでいるときにのみ、行$I$, 列$J$が1になる。
なお、これを合成行列と呼び、こういうモデルを合成リンクモデルという。提唱者はThompson & Baker (1981)で、一般化線形モデルの拡張になっている。[←そうなの？ よくわかんないけど、まあいいや]
$\mathbf{\gamma}$は非負じゃないと困るので、
$\mathbf{\gamma} = \exp(\mathbf{X} \mathbf{\beta})$
とし、$\mathbf{\beta}$ を推定する。普通は、$\mathbf{X}$は$J \times J$の単位行列、$\mathbf{\beta}$は長さ$J$のベクトル、でよい。以下でもそう考えるが、$J$がすっごく大きいときはパラメータ数が大きくなりすぎるので、$\mathbf{X}$を次元$p$のBスプライン基底としてもよい。
さて、こうしてデータ生成プロセスを妄想するのは勝手だが、問題は、$\mathbf{\beta}$をどうやって推定するかである。
$\mathbf{\gamma}$はスムーズだと考え、$\mathbf{\beta}$の非スムーズさを表す行列$\mathbf{D}_2$を求める[詳細は付録。読んでないんだけど、ちらっとみたところでは、ある$\beta_j$について、その両隣の和からその2倍を引いた値を求める模様]。で、ペナルティ
$P = (\mathbf{D}_2 \mathbf{\beta})^2$
というのを考える。
さて、このモデルのポアソン対数尤度は
$l = \sum_i^I (y_i \log\mu_i -\mu_i)$
である。ここからペナルティを引いた罰則つき対数尤度
$l^* = l - \frac{\lambda}{2} P$
を最大化する。これはiteratively reweighted least-squaresで解ける。
パラメータ$\lambda$は、AICが最小となる値をグリッドサーチして決めるがよかろう。
実データへの適用例... [省略]
考察。
本提案は頻度主義アプローチである。前にこれのベイジアン版を提案したが(そのときはオープンエンドのことを考えてなかった)、それだと$\lambda$の決定に伴う不確実性を考慮できる。シミュレーションで比較してみたが(付録)、たいした違いはなかった。
このモデルの漸近特性は解析的にはわかっておらず、今後の課題。シミュレーションで調べたところ(付録)、一致性など良い性質がみられた。
今後の課題は...まず2次元への拡張。各ビンにおける平均やSDが分かっているときにそれをうまく使う方法。
... なるほどねー。データの背後に平滑な離散分布を考えるわけね。そこでのカテゴリ数($J$)はどのくらいがよいのだろうか？ あまり大きくするとなんらかの構造($\mathbf{X}$)を入れないといけなくなるだろうけど、小さすぎても良くないような気がする。
論文：データ解析(2018-) - 読了: Rizzi, Gampe, & Eilers (2015) ヒストグラムから確率分布をノンパラメトリック推定する罰則つき合成リンクモデル