« 読了：Reyes, Francisco-Fernandez, & Cao (2016) ヒストグラムから元の確率密度を推定します、階級の幅が不均等でも大丈夫です | メイン | 読了: Rizzi, Gampe, & Eilers (2015) ヒストグラムから確率分布をノンパラメトリック推定する罰則つき合成リンクモデル »
2018年10月16日 (火)
ヒストグラムから確率密度曲線をノンパラメトリックに復元し、任意の区間の確率を推定したい... 絶対そういうのがあるはずだ... と、ディスプレイの前で延々と、探し、探し求めてー [←突然「長崎は今日も雨だった」の節になる]、ついに見つけた論文。
Rizzi, S., Thinggraard, M., Engholm, G., Christensen, N., Johannessen, T.B., Vaupel, J.W., Lindahl-Jacobsen, R. (2016) Comparison of non-parametric methods for ungrouping coarsely aggregated data. BMC Medical Research Methodology. 16:59.
いやー、正直いって、この論文に辿り着くまでに半日かかった。あきらかに探し方が悪かったんだけど、ポジティブにいえば、その半日で少しは賢くなったということだろう。
著者らの整理によれば、その方法は大きく３つある。
その1, カーネル密度推定量。
いま、未知の連続分布$f(x)$からの標本$x_1, \ldots, x_n$があるとして、$f(x)$のカーネル密度推定量は
$\hat{f}(x) = \frac{1}{nh} \sum_i^n K(\frac{x-X_i}{h})$
この推定においてポイントになるのは、カーネル関数$K()$よりもむしろバンド幅$h$である。$h$はプラグイン推定かCVで決めることが多い。
さて、ヒストグラムのように値がグループ化されてて、もはや$x_1, \ldots, x_n$がわからないときはどうするか。
- Scott & Sheather (1985): 初期の提案。すべて階級の中央にあったことにしちゃう。階級幅が細かいときにはこれでもうまくいく。
- Blower & Kelsall(2002): $K(x-x_{ji})$($j$は階級番号)のかわりに$E[K(x-x_{ji})]$として、反復計算でなんとかする [←よくわからん]
- Wang & Wertelecki (2013): ブートストラップ型カーネル密度推定量を提案。
- Braun et al.(2005) 条件付き期待値に基づく局所尤度推定。
その2, スプライン補完。
ヒストグラムを累積分布に直し、データ点の間を2次だか3次だかのスプライン関数でつなぐ。よく用いられている方法なんだけど、オープンエンドの階級のときに困る。うっかりすると曲線が下がっちゃう(つまりある区間の頻度が負になる)から。対策:
- Human Mortality Databaseというのがあって、そこではオープンエンドのところだけパラメトリックモデルを使っている。
- Smith et al.(2004): 3次スプライン補完に単調制約をかける。Hymanフィルタという由。
- Human Fertility Databaseというのもあって、そこではpiecewise 3次Hermite補完多項式なるものを使っておる。
その3, penalized composite linkモデル。
composite linkモデルというのは一般化線型モデルの拡張で、$I$個の階級について頻度$y_i$がわかっているとき、それを期待値$\mu_i$のポアソン分布の実現値とみる。で、実は背後にもっと狭い$J$個の階級があり、それぞれが$\gamma_j$を持っていて、それを$I$個の階級にまとめたら$\mu_i$が得られたのだと考える。この$\gamma_j$の分布を最尤推定するんだけど、その際に曲線のラフさを表す罰則項をいれる(そうしないと解けない)。この方法でポイントになるのは罰則項に掛けるパラメータで、AICで決める。
さて、文献を検索しまくりまして、この問題に使われてきた方法を探しました。次の5つが見つかった。
- bootkde:ブートストラップ・カーネル密度推定量。Rでいうと、bda::bde()の"bootkde"メソッド。
- hermite spline: ピースワイズ3次Hermite補完多項式。signal::interp1()の"pchip"メソッド。
- hyman spline: Hymanフィルタ付きスプライン補完。demography::cm.spline()。
- ickde: iterated conditional expectationカーネル密度推定量。ICD::ickde()。バンド幅を決める方法として、プラグイン推定はKernSmooth::dpik()、CVはICE::bickde()がある。
- pclm: penalized composite linkモデル。Rizzi et al. (2015 Am.J.Epidemiol.)がコードを公開している。
オーケー、じゃシミュレーションしてみましょう。
架空の年齢別死亡数の分布を考える。あるパラメータ[略]を持つワイブル分布だとする。ここから標本を抽出し、ヒストグラムをつくり、真の分布を推定する。これを500回繰り返す。
動かす条件は、標本サイズ{200, 1000}, 階級の切り方{5年幅, 5年幅だが85歳より右はオープンエンド}。2x2=4つのシナリオができるわけだ。
成績の指標は3つ。(1)$IAE = \sum|\hat{f}(x)-f(x)|$。(2)$ISE = \sum(\hat{f}(x)-f(x))^2$。(3)KL距離 $\sum f(x) \log(f(x)/\hat{f}(x))$。
他にNORDCANデータベースというのから、実際の死亡データを貰ってきて分析してみた[このくだりは略]。
結果。
ickdeでバンド幅を決める方法はKernSmooth::dpik()がよかったので、以下そっちを使う。
オープンエンドなしの場合だと、n=200ではスプライン補完(hermite spline, hyman spline)の成績が悪い。n=1000では大差なし。
オープンエンドがあると、bootkde, hermite spline, ickdeでは、推定された密度曲線が右の方で変な風に跳ねる[←へぇー]。pclmが最優秀, hyman splineがこれに次ぐ。
[実データの話。略]
考察。
penalized composite linkモデルが最優秀、Hymanフィルタ付きスプライン補完がこれに次ぐ、という結果になりました。
... ヒストグラムにオープンエンドの階級があると、カーネル密度推定が上手くいかないって話、面白いなあ。素人ながら想像するに、これはカーネル密度推定の本質的な問題というより、オープンエンドをどう扱うかという問題であって、改善の余地がある話なのかもしれない。
いま気がついたんだけど、この論文の第一著者のRizziさん、penalized composite linkモデルのコードを公開している人なのね。自分の提案手法を推していたのか、ははは。調べてみたところ、penalized composite linkモデルのRパッケージとしてungroupというのがあって、開発者のひとりがRizziさん。この論文のあとで公開したのだろう。
検討されている中でいちばんシンプルだなと思うのは、累積分布を単調制約付きでスプライン補完するというアイデアである。標準パッケージstatsにsplinefun()というのがあって、hyman法というのを指定できるんだけど、これじゃだめなんだろうか。
ところで... ヒストグラムがセンサスに由来している場合、推定された累積密度曲線は、ヒストグラムの階級境界においては、もとの累積分布と交わらないと気持ち悪いような気がする。こういう場合はスプライン補完一択、ってことなのかな？
論文：データ解析(2018-) - 読了：Rizzi et al. (2016) ヒストグラムから元の確率分布をノンパラメトリックに推定する方法コンテスト