kur0cky

こんにちは

「大規模計算時代の統計推論」 を全部Rでやってみる ~第4章~

こんにちは,kur0ckyです. 夏もやっと終わりが見え,涼しい日には好きな長袖が着れて嬉しいです. 秋はめっちゃ好きです.

さて,今回扱う第4章は「フィッシャー派的な推論と最尤推定」です. 3章までは,頻度派とベイズ派の二軸で説明されてきましたが,ここで新たに「フィッシャー派」というもう一つの流派が出てきました. また,4章のはじめには次のような記述があります.

「分散分析」「有意差検定」「最尤推定」は,ほとんどいつも頻度派的に適用されてきた.しかし,それらフィッシャー派の手法の理論的根拠の多くは,実際にはベイズ派と頻度派のいずれの発想にもよらないものであり,時には両者の組み合わせによるものであった.

しっかり読んでいくことで,「フィッシャー派」がどのような立ち位置にあるのか,理解したいです. できないかもしれません. それではやっていきましょう.

4.1 尤度と最大尤度

まずは最尤推定のおさらいです. 基本的な概念をさらった後,腎臓病患者のgfrデータを元にいくつかの分布のパラメータを最尤推定する例が紹介されています.

基本的な概念

なんらかの確率モデルを仮定したとき,そのモデルのもと観測値 \bf xが得られる確率を尤度 (likelihood)といいます. また,観測値を定数とし,パラメータ \muの関数と捉えたものを尤度関数 (likelihood function)といったりします. つまり,尤度関数 L _ {\bf x}(\mu)は次で表されます.

 L _ {\bf x}(\mu) = p({\bf x} | \mu).

さらに, \bf xの各要素が独立に観測されるとしたとき,尤度関数は単純な積の形に分解できます.

 L _ {\bf x}(\mu) = p({\bf x} | \mu) = \prod _ {i} p(x _ i | \mu).

この尤度関数の自然対数をとったものを,対数尤度関数(log likelihood function)といい,ここでは l _ {\bf x} = \log L _ {\bf x}(\mu)と表します. 対数をとることにより,和の形で表すことができます. 積の形のままでは, L _ {\bf x}(\mu)があっという間に小さくなってしまうことがありますが,和の形ではこれを防ぐことができます. もちろん,単純に扱いやすいこともメリットです. そして,この対数尤度を最大にするようなパラメータの値が最尤推定値 (maximum likelihood estimate; MLE)です. モデルのパラメータの関数になっているような量に興味がある場合は,もちろん最尤推定値を差し込んでも良いです.

 \text{MLE : }\hat\mu = \text{argmax} _ {\mu \in \Omega} l _ {\bf x}(\mu).

ちなみに,前回の記事でも紹介しましたが,最尤推定値は無情報事前分布を用いた場合のMAP推定値(事後分布を最大にするような推定値)と同じになります. 無情報事前分布は一様(定数)であるためです.

 
\mu _ {MAP} =
\text{argmax} _ {\mu} p(\mu|{\bf x}) =
\text{argmax} _ {\mu}p({\bf x}|\mu)p(\mu) = 
\text{argmax} _ {\mu}p({\bf x} | \mu)

そのため,本書でも

最尤推定は妥当なベイズ派の正統性も有している.

という記述があります.

図4.1(正規分布とガンマ分布の最尤推定

さっそく解析例です. 2章でも扱った腎臓病患者のgfrデータをもとに,正規分布とガンマ分布の母数を推定してみましょう. 通常ガンマ分布の下限は x \geq 0なのですが,本書では様々な下限 \lambdaをとりうるようにした分布に対して推定を行っています. この特殊なガンマ分布は次式で表されます.

 
\displaystyle
f _ {\mu} (x) = 
\begin{cases}
\frac{(x - \lambda) ^ {(x - \nu)}}{\lambda ^ {\nu}\Gamma(\nu)}\exp(-\frac{x-\lambda}{\sigma}) & x > \lambda, \\\
0 & x \leq \lambda.
\end{cases}

正規分布やガンマ分布の場合は,最尤推定値が解析的に求まるのですが,この特殊な場合は求まらないそうです. この記事では,面倒なので全部数値的に求めてしまいます.笑

それではやっていきましょう. まずはパッケージを読み込み,データをダウンロードします.

library(tidyverse)
library(stats4)

theme_set(theme_classic(base_family = "HiraKakuPro-W3"))

gfr <- scan("https://web.stanford.edu/~hastie/CASI_files/DATA/gfr.txt")

Rには最尤推定をするための様々なパッケージがありますが,ここではstats4パッケージのmle()を使ってみることにします. stats4は追加のインストールが必要なく,歴史も長いです. アルゴリズムはシンプルに準ニュートン法を使うことにします.

mle()には,負の対数尤度関数を引数に要求し,これを最小化します(つまり対数尤度の最大化です). 対数尤度関数を用意しましょう. せっかくなので普通のガンマ分布も推定します. 正規分布と普通のガンマ分布はそれぞれdnorm(), dgamma()が密度関数として用意してあるので非常に簡単です. 下限が0でないガンマ分布の対数尤度は,少し計算して次のようになります(ただし x > \lambdaで).

 \displaystyle
l _ {\bf x}(\lambda, \sigma, \mu) = 
\sum _ i
\left (
 (\nu -1)\log (x _ i - \lambda) - \nu \log(\sigma) - \log(\Gamma(\nu)) - (x _ i -\lambda) / \sigma
\right) .
# 正規分布
neg_LL_normal <- function(mu,sigma) {
  x = gfr
  -sum(dnorm(x, mean=mu, sd=sigma, log = T)) 
}

# 普通のガンマ分布
neg_LL_gamma <- function(shape, rate){
  x = gfr
  -sum(dgamma(x, shape, rate,  log = T))
}

# 下限が自由なガンマ分布
neg_LL_gamma2 <- function(lambda, sigma, nu){
  x  = gfr
  ll = (nu - 1) * log(x - lambda) - nu * log(sigma) - log(gamma(nu)) - (x - lambda) / sigma
  ll[x < lambda] = 0
  -sum(ll)
}

準備ができたので最適化していきましょう. 数値的な発散によるエラーを避けるため,少しズルをして初期値をうまく与えています. 本当なら,様々な初期値を用意して,最も対数尤度が大きくなる値を選ぶべきでしょう.

normal_fit <- mle(minuslogl = neg_LL_normal,
                  start=list(mu=50, sigma=10),
                  method = "L-BFGS-B",
                  upper = c(Inf, Inf),
                  lower = c(-Inf, 0))
gamma_fit <- mle(minuslogl = neg_LL_gamma,
                 start=list(shape = 50, rate = 0.1),
                 method = "L-BFGS-B")
gamma2_fit <- mle(minuslogl = neg_LL_gamma2,
                 start=list(lambda = 21.5, sigma = 5.5, nu = 6),
                 method = "L-BFGS-B")

最後に図を書いて終わらせます. 推定値はfit@coefに入っています.

colors <- c(normal = "black",
            gamma = "blue",
            truncated_gamma = "steelblue")
tibble(gfr = gfr) %>% 
  ggplot()+
  geom_histogram(aes(gfr, y = ..density..), 
                 colour = "black", fill = "green", size = .3)+
  stat_function(fun = dgamma2, aes(colour = "truncated_gamma"),
                size = 1, linetype = 2,
                args = list(lambda = gamma2_fit@coef["lambda"],
                            sigma = gamma2_fit@coef["sigma"],
                            nu = gamma2_fit@coef["nu"]))+
  stat_function(fun = dgamma, aes(colour = "gamma"),
                size = 0.5,
                args = list(shape = gamma_fit@coef["shape"],
                            rate = gamma_fit@coef["rate"]))+
  stat_function(fun = dnorm, aes(colour = "normal"), size = 1,
                args = list(mean = normal_fit@coef["mu"],
                            sd = normal_fit@coef["sigma"]))+
  scale_colour_manual(values = colors)+
  scale_x_continuous(limits = c(15, 118))+
  labs(x = "gfr", y = "頻度")

綺麗な図がかけました. 今回は関数のプロットにstat_function()を使いました. stat_function()は,引数argsにリストを渡す形でパラメータを指定することができます. ちなみに,この中では尤度,AICともに下限が自由なガンマ分布が良くなるようです.

f:id:kur0cky:20200921162504p:plain

4.2 フィッシャー情報量と最尤推定

フィッシャー派の大本命であろうと思われるフィッシャー情報量について書かれた節です. 解析がないのですっ飛ばします.

フィッシャー情報量については前回の記事でもチラッとだけ紹介しました(ほぼ他の記事に丸投げでしたが).

kur0cky.hatenablog.com

4.3 条件付き推論

この節では条件付き推論 (conditional inference) と補助統計量 (ancillary statistic) について書かれています. 補助統計量とは,着目するパラメータ \muと無関係な分布をもつ確率変数や統計量のことです. たとえば,この節の冒頭で出てくるコイン投げの確率変数などです.  \muに関する推論を行う際,このような確率変数はある値に固定してしまってよいことになります.

この記事で扱うには少し重すぎるので,詳しくはwikipediaなどをご覧ください. 余裕があれば,十分統計量やBasuの理論などと一緒に解説する記事を書きたいと思います.笑

en.wikipedia.org

図4.2(コーシー分布のフィッシャー情報量の例)

補助統計量の例のうち,コーシー分布の例は解析があるので,この記事でも再現します. この例を通して,観測されたデータから求めたフィッシャー情報量が近似的な補助統計量として振舞うことが感覚的に示されています. つまり, コーシー分布のフィッシャー情報量を計算しましょう.

解析の手順は大きく次の5ステップです

  1. コーシー分布から大きさ20の標本を10000個発生させる.
  2. それぞれの標本に対し,locationの最尤推定 \hat \thetaを求める.
  3. また,それぞれの標本に対し,観測情報量下限 1/I({\bf x})を計算する.
  4.  1/I({\bf x})の10分位点ごとに \hat\thetaの経験分散を求める.
  5. 分位点ごとの経験分散と,分位点ごとの中央値をプロットする.

ここで, I({\bf x})最尤推定量を差し込んだフィッシャー情報量を使った n\mathcal I _ {\hat\theta}ではなく,最尤推定値を代入したもので,次式で表されます.

 \displaystyle
I({\bf x}) = - \ddot l _ {\bf x} (\hat \theta) = - \frac{\partial ^ 2}{\partial\theta ^ 2} l _ {\bf x} (\theta) | _ {\hat\theta}.

まずは,標本から「観測されたフィッシャー情報量」を計算する関数を作ります. まずは対数尤度の二回微分をしなければなりません. ゴリゴリ数式で攻めても良いのですが,全てRで解析する企画なので,微分もRでやってしまいます(楽してるわけではありません笑). Rでは,expression()で対数尤度関数の数式オブジェクトを作っておき,D()微分することができます.

f <- expression(log(1/pi/(1 + (x-theta)^2)))

D(D(f, "theta"), "theta")

綺麗に整理するところまではやってくれませんが,微分された数式が出てきます. これをコピペし,dd_ll_cauchy()という関数を作っておきましょう.

dd_ll_cauchy <- function(x, theta){
  -((1/pi * 2/(1 + (x - theta)^2)^2 - 1/pi * (2 * (x - theta)) * (2 * (2 * (x - theta) * (1 + (x - theta)^2)))/((1 + (x - theta)^2)^2)^2)/(1/pi/(1 + (x - theta)^2)) + 1/pi * (2 * (x - theta))/(1 + (x - theta)^2)^2 * (1/pi * (2 * (x - theta))/(1 + (x - theta)^2)^2)/(1/pi/(1 + (x - theta)^2))^2)
}

続いて標本を発生させ,location  \theta最尤推定,そして「観測されたフィッシャー情報量」の算出まで一気にやってしまいます.

set.seed(1)
df <- tibble(index = 1:10000) %>% 
  mutate(x = map(index, ~  rcauchy(20, 0, 1))) %>% 
  mutate(theta_hat = map_dbl(x, ~ 
                               mle(minuslogl = function(theta){-sum(dcauchy(.x, location = theta, scale = 1, log = T))},
                                   start=list(theta = 0),
                                   method = "L-BFGS-B") %>% 
                               .@coef)) %>% 
  mutate(f_info = map2_dbl(x, theta_hat, ~ 
                           -1/sum(dd_ll_cauchy(.x, .y))))
print(df)
# A tibble: 10,000 x 4
   index x          theta_hat f_info
   <int> <list>         <dbl>  <dbl>
 1     1 <dbl [20]>   -0.132  0.113 
 2     2 <dbl [20]>    0.120  0.128 
 3     3 <dbl [20]>   -0.512  0.126 
 4     4 <dbl [20]>    0.0721 0.140 
 5     5 <dbl [20]>   -0.0773 0.183 
 6     6 <dbl [20]>    0.111  0.0983
 7     7 <dbl [20]>   -0.0705 0.108 
 8     8 <dbl [20]>    0.413  0.108 
 9     9 <dbl [20]>   -0.384  0.0799
10    10 <dbl [20]>   -0.118  0.0960
# … with 9,990 more rows

求めたフィッシャー情報量を100分率の区間 (0-5), (5-15), \dots, (95-100)の分位点にわけ, \hat\thetaの経験分散を求めましょう. percent_rank()を使うとパーセンタイルが求まるので,それを100倍したものが区切りbreaksを上回っている数でグルーピングします.

# 10分位点の区切り
breaks <- c(seq(5, 95, 10), Inf)

df <- df %>% 
  mutate(percentile = percent_rank(f_info) * 100,
         rank = map_int(percentile, ~ sum(.x > breaks))) %>% 
  group_by(rank) %>% 
  summarise(f_info = mean(f_info),
            se = var(theta_hat)) %>%
  drop_na() 
print(df)
# A tibble: 11 x 3
    rank f_info     se
   <int>  <dbl>  <dbl>
 1     0 0.0547 0.0614
 2     1 0.0648 0.0682
 3     2 0.0728 0.0782
 4     3 0.0794 0.0830
 5     4 0.0863 0.0969
 6     5 0.0937 0.0997
 7     6 0.102  0.107 
 8     7 0.113  0.111 
 9     8 0.128  0.134 
10     9 0.157  0.186 
11    10 0.273  0.242 

最後に作図です. 今回は単純な散布図なのでコードも少なく,楽でした.

df %>% 
  ggplot(aes(f_info, se))+
  geom_point()+
  geom_abline(slope = 1, linetype = 2)+
  geom_hline(yintercept = 0.1, linetype = 2, colour = "red")+
  labs(x = "観測された情報量下限", y = "最尤推定量の分散")

f:id:kur0cky:20200915213030p:plain

疑問点

最尤推定量は一致推定量ではありますが,不偏推定量ではないですよね?? クラメル=ラオの下限は不偏推定量の中での話なので,この例のような比較を行うことには意味があるのか疑問に思いました. 大きさ20のような小標本でも,今回の例のように,下限に近い値はでますよ,ということが言いたかったのでしょうか. 大標本なら漸近的に不偏になるので良さそうですが.

おまけ:数式でやる

どうせなので,コーシー分布のlocationのフィッシャー情報量を数式でみましょう. コーシー分布の尤度関数は次式で与えられます.

 \displaystyle
f({\bf x}| \theta) = \prod\frac{1}{\pi(1 + (x _ {i} -
\theta) ^ 2)}.

対数尤度を \theta微分したスコア関数は,素直に計算し次のようになります.

 \displaystyle
\frac{\partial}{\partial \theta} \log f({\bf x}| \theta) = 
\sum \frac{\frac{2}{\pi(1 + (x _ {i} -
\theta) ^ 2) ^ 2}(x-\theta)}{f(x _ {i}| \theta)} = 
\sum \frac{2(x _ i - \theta)}{1 + (x _ i - \theta)^ 2}.

次に,スコア関数の二次モーメントは次のように計算できます. 途中で u = x-\theta, v = 1/(1 + u ^ 2)の置換積分をしています. 最後はベータ関数の形を利用して整理しています.

 \displaystyle
\begin{align}
\mathcal I _ {x _ i} (\theta)  &= 
{\rm E}\left(\left(\frac{\partial}{\partial \theta} \log f(x _ i| \theta)\right) ^ 2\right) \\\ &= 
\int _ {-\infty} ^ {\infty} \left(\frac{2(x _ i - \theta)}{1 + (x _ i - \theta)^ 2}\right) ^ 2 \frac{1}{\pi(1 + (x-\theta)^ 2)}dx \\\ &=
\frac{4}{\pi} \int _ {-\infty} ^ {\infty} \frac{(x-\theta) ^ 2}{(1 + (x-\theta) ^ 2) ^ 3}dx \\\ &=
\frac{8}{\pi}\int _ {0} ^ {\infty}\frac{u ^ 2}{(1 + u ^ 2) ^ 3}du \\\ &= 
\frac{4}{\pi}\int _ 0 ^ 1\sqrt{v}\sqrt{1-v} dv \\\ &=
\frac{4}{\pi}\int _ 0 ^ 1 v ^ {3/2 - 1} (1-v) ^ {3/2 - 1} dv \\\ &=
\frac{4\Gamma(3/2)\Gamma(3/2)}{\pi\Gamma(3/2 + 3/2)} \\\ &=
\frac{4(0.5\sqrt{\pi}) ^ 2}{\pi 2!} = \frac{1}{2}.
\end{align}

よって,フィッシャー情報量はこの合計であり,

\displaystyle
 \mathcal I _ {\bf x} (\theta) = \frac{n}{2}.

となります.ありがとうございました. 前節での例では,大きさ20の標本を作っているので,クラメルラオの下限は0.1となります.

4.4 並び替えと無作為化

 t検定をはじめとする多くの手法は,標本が正規分布に従うことを仮定します. この仮定はしばしば現実てきでなく,批判にさらされます. 本書では,これに対するフィッシャーの解決策として,並べ替えを利用したノンパラメトリック検定が紹介されています.

図4.3

まずはデータセットをダウンロードします. ついでに,対象の遺伝子である136行目を抜き出します.

leukemia <- read_csv("http://web.stanford.edu/~hastie/CASI_files/DATA/leukemia_big.csv")

leukemia_136 <- leukemia[136, ] %>% 
  unlist() %>% 
  unname()

シャッフルしないもとの t検定は1章でやっていたので,ここでは省略します.

続いて,データをシャッフルし47個の値と25個の値に分け,これに対して t検定します. これを B = 10000回やります. map()nest()が大好きなので次のように実装しました.

B <- 10000
set.seed(1)

df_sim <- tibble(index = 1:B) %>% 
  mutate(split_index = map(index,~ sample(1:72, 25, replace = FALSE)),
         x1 = map(split_index, ~ leukemia_136[.x]),
         x2 = map(split_index, ~ leukemia_136[-.x]),
         t = map2_dbl(x1, x2, ~ t.test(.x, .y)$statistic))

print(df_sim)

次のようなtibble()ができます. 何度でも推したいですがnest()は神です.

# A tibble: 10,000 x 5
   index split_index x1         x2              t
   <int> <list>      <list>     <list>      <dbl>
 1     1 <int [25]>  <dbl [25]> <dbl [47]> -1.66 
 2     2 <int [25]>  <dbl [25]> <dbl [47]>  0.662
 3     3 <int [25]>  <dbl [25]> <dbl [47]>  3.07 
 4     4 <int [25]>  <dbl [25]> <dbl [47]>  0.276
 5     5 <int [25]>  <dbl [25]> <dbl [47]> -0.950
 6     6 <int [25]>  <dbl [25]> <dbl [47]> -0.411
 7     7 <int [25]>  <dbl [25]> <dbl [47]>  0.752
 8     8 <int [25]>  <dbl [25]> <dbl [47]>  0.381
 9     9 <int [25]>  <dbl [25]> <dbl [47]> -0.853
10    10 <int [25]>  <dbl [25]> <dbl [47]>  0.895
# … with 9,990 more rows

ALLAMLの間に差がないとする帰無仮説に対して有意水準を求めます. sum(abs(df_sim$t) >= 3.01) / Bで求めると0.0045となり,本書の0.0025とは少し離れてしまいましたが,乱数のせいにしておきましょう.

続いて作図します. annotate()で注釈や矢印を,geom_segment()でx軸のポイントを ちまちまと書いていきます. 軸ラベルのカスタマイズは,scale_x_continuous()内のbreaksで指定します.

df_sim %>% 
  ggplot(aes(t))+
  geom_histogram(fill = "gray100", colour = "steelblue",
                 bins = 50)+
  geom_vline(xintercept = -3.01, linetype = 2, colour = "red")+
  geom_vline(xintercept = 3.01, linetype = 2, colour = "red")+
  annotate("text", x = 3.5, y = 100, label = "もとのt統計量", family = "HiraKakuPro-W3")+
  annotate("segment", x = 3.5, xend = 3.1, y = 70, yend = 20,
           size=0.5, arrow=arrow(length=unit(0.30,"cm")))+
  scale_x_continuous(breaks=c(-4, -3.01, -2, 0, 2, 3.01, 4))+
  geom_segment(data = filter(sim, abs(t) > 3.01),
               aes(x = t, xend = t, y = 0, yend = -20),
               colour = "red", size = 0.3)+
  labs(x = "t*値",
       y = "頻度")

f:id:kur0cky:20200915213034p:plain

このようなシミュレーションによる検定は,正規分布の仮定をしなくても, X _ i \sim f _ \mu (x)の独立同分布を帰無仮説として行うことができます. いわゆるノンパラメトリック検定の一つですね. ちなみに,この実験からは,ALLの遺伝子136の値に比べAMLの値が大きいと言えるかを判断することができますが,因果までは言うことができません. これと対照的な,有名なフィッシャーの農場実験の話題にも本書では触れられています.

疑問点

ただ,一つだけ疑問を感じたので,書いておきます. たしかに,帰無仮説 X _ i \sim f _ \mu (x)のもとでは,シャッフルして分割した 72!/(42!25!)通りの組み合わせは同様に確からしいです. しかし, f _ \mu (x)を経験分布で近似すると考えると,並び替えではなく無作為に復元抽出するべきではないでしょうか. どちらが良いか私の弱い知識では判断できなかったので,どなたか教えてください(懇願)...

おわりに

内容もだんだんと骨太になってきました. 少ししんどいですが細々とやっていきます. ではでは