進化学実習 2022 牧野研

6. 統計モデリング1: 確率分布、尤度

岩嵜航
東北大学 生命科学研究科 進化ゲノミクス分野 牧野研
2022-04-11 東北大学 理学部生物学科 進化学実習 https://heavywatal.github.io/slides/tohoku2022r/

データを使ってやりたいこと

  • 現象を理解したい
  • 将来を予測したい
  • ものを分類・判別したい
  • 挙動を制御したい
  • 新しい何かを生成したい

そのために解析は必要? 未加工の生データこそ宝?

データ解析って必要? 生データ見ればいいべ?

往々にして複雑過ぎ、情報多すぎ、そのままでは手に負えない

print(ggplot2::diamonds)
      carat       cut color clarity depth table price     x     y     z
      <dbl>     <ord> <ord>   <ord> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
    1  0.23     Ideal     E     SI2  61.5    55   326  3.95  3.98  2.43
    2  0.21   Premium     E     SI1  59.8    61   326  3.89  3.84  2.31
    3  0.23      Good     E     VS1  56.9    65   327  4.05  4.07  2.31
    4  0.29   Premium     I     VS2  62.4    58   334  4.20  4.23  2.63
   --                                                                  
53937  0.72      Good     D     SI1  63.1    55  2757  5.69  5.75  3.61
53938  0.70 Very Good     D     SI1  62.8    60  2757  5.66  5.68  3.56
53939  0.86   Premium     H     SI2  61.0    58  2757  6.15  6.12  3.74
53940  0.75     Ideal     D     SI2  62.2    55  2757  5.83  5.87  3.64

ダイヤモンド53,940個について10項目の値を持つデータセット

要約統計量を見てみよう

各列の平均とか標準偏差とか:

   stat carat depth table    price     x     y     z
  <chr> <dbl> <dbl> <dbl>    <dbl> <dbl> <dbl> <dbl>
1  mean  0.80 61.75 57.46  3932.80  5.73  5.73  3.54
2    sd  0.47  1.43  2.23  3989.44  1.12  1.14  0.71
3   max  5.01 79.00 95.00 18823.00 10.74 58.90 31.80

大きさ carat と価格 price相関係数は 0.92 (かなり高い)。

生のままよりは把握しやすいかも。

しかし要注意…

平均値ばかり見て可視化を怠ると構造を見逃す

https://www.autodesk.com/research/publications/same-stats-different-graphs

作図してみると全体像・構造が見やすい

情報をうまく絞って整理 → 直感的にわかる

plot of chunk simplify-diamonds

carat が大きいほど price も高いらしい。
その度合いは clarity によって異なるらしい。

統計とは

データをうまくまとめ、それに基づいて推論するための手法。

  • 記述統計: データそのものを要約する
    • 要約統計量 (e.g., 平均、標準偏差、etc.)
    • 作図、作表
  • 推測統計: データの背後にある母集団・生成過程を考える
    • 数理モデル
    • 確率分布
    • パラメータ(母数)

「グラフを眺めてなんとなく分かる」以上の分析にはモデルが必要

モデルとは

対象システムを単純化・理想化して扱いやすくしたもの

Mathematical Model 数理モデル
数学的な方程式として記述されるもの。
e.g., Lotka-Volterra eq., Hill eq.

Computational Model 数値計算モデル
数値計算の手続きとして記述されるもの。
e.g., Schelling’s Segregation Model, tumopp

Concrete Model 具象モデル
具体的な事物で作られるもの。
e.g., San Francisco Bay-Delta Model
Weisberg 2012 "Simulation and Similarity" (科学とモデル)

データ科学における数理モデル

データ生成をうまく真似できそうな仮定の数式表現。
 


「データ分析のための数理モデル入門」江崎貴裕 2020 より改変

データ科学における数理モデル

データ生成をうまく真似できそうな仮定の数式表現。
e.g., 大きいほど高く売れる: $\text{price} = A \times \text{carat} + B + \epsilon$

plot of chunk lm-diamonds

新しく採れたダイヤモンドの価格予想とかにも使える。

このように「YをXの関数として表す」ようなモデルを回帰と呼ぶ。

今回は回帰を軸とした統計モデリングの解説

単純な直線あてはめから出発し、ちょっとずつ統計モデリング。

久保さん https://kuboweb.github.io/-kubo/ce/LinksGlm.html

回帰は教師あり機械学習の一種とも言える

でも統計モデリングはいわゆる“機械学習”とは違う気もする…?

モデリングにおける2つのアプローチ


「データ分析のための数理モデル入門」江崎貴裕 2020 より改変

どっちも知っておいて使い分けたい

項目 統計モデリング 近年の機械学習
一般化線形モデル
階層ベイズモデル
ランダムフォレスト
ニューラルネットワーク
モデル構造 単純化したい 性能のためなら複雑化
モデル解釈 ここが強み 難しい。重視しない。途上。
予測・生成 うまくすれば頑健 主目的。強力。高精度
データ量 少なくてもそれなり 大量に必要
計算量 場合による 場合による

統計モデリングのお品書き

データ解析のための統計モデリング入門 久保拓弥 2012

久保先生の"緑本"こと
「データ解析のための統計モデリング入門」
をベースに回帰分析の概要を紹介。

  1. イントロ (#6 今回)
  2. 統計モデルの基本
    • 確率変数・確率分布 👈 主役
    • 尤度・最尤推定
  3. 一般化線形モデル、混合モデル (#7 次回)
  4. ベイズ統計、階層ベイズモデル 今回はここまでやらない

回帰のキモは線ではなく分布


Data Science Hill Climb 2021 (東京海上) での講義 (~6時間) の演習無し抜粋バージョン (~2時間 x 2回)。

回帰モデルの2段階

  1. Define a family of models: だいたいどんな形か、式をたてる

    • 直線: $y = a_1 + a_2 x$
    • 対数: $\log(y) = a_1 + a_2 x$
    • 二次曲線: $y = a_1 + a_2 x^2$
  2. Generate a fitted model: データに合うようにパラメータを調整

    • $y = 3x + 7$
    • $y = 9x^2$
https://r4ds.had.co.nz/model-basics.html

たぶん身長が高いほど体重も重い

なんとなく $y = a x + b$ でいい線が引けそう
 

plot of chunk weight-height

たぶん身長が高いほど体重も重い

なんとなく $y = a x + b$ でいい線が引けそう
じゃあ切片と傾き、どう決める?

plot of chunk weight-lines

最小二乗法

回帰直線からの残差平方和(RSS)を最小化する。

plot of chunk weight-residual

残差平方和(RSS)が最小となるパラメータを探せ

ランダムに試してみて、上位のものを採用

plot of chunk weight-goodlines

残差平方和(RSS)が最小となるパラメータを探せ

グリッドサーチ: パラメータ空間の一定範囲内を均等に試す

plot of chunk weight-grid

こうした最適化の手法はいろいろあるけど、ここでは扱わない。

これくらいなら一瞬で計算してもらえる

par_init = c(intercept = 0, slope = 0)
result = optim(par_init, fn = rss_weight, data = df_weight)
result$par
intercept     slope 
-50.54532  67.18659 

plot of chunk weight-lm

上記コードは最適化一般の書き方。
回帰が目的なら次ページのようにするのが楽 →

lm() で直線あてはめしてみる

fit = lm(data = mpg, formula = hwy ~ displ)
broom::tidy(fit)
         term  estimate std.error statistic       p.value
        <chr>     <dbl>     <dbl>     <dbl>         <dbl>
1 (Intercept) 35.697651 0.7203676  49.55477 2.123519e-125
2       displ -3.530589 0.1945137 -18.15085  2.038974e-46
mpg_added = modelr::add_predictions(mpg, fit)
ggplot(mpg_added) + aes(displ, hwy) + geom_point() +
  geom_line(aes(y = pred), size = 1, color = "#3366ff")

plot of chunk lm-mpg

🔰 diamonds などほかのデータでも lm() を試してみよう。

何でもかんでも直線あてはめではよろしくない

plot of chunk lm-bad

  • 観察データは常に正の値なのに予測が負に突入してない?
  • 縦軸は整数。しかものばらつきが横軸に応じて変化?

何でもかんでも直線あてはめではよろしくない

plot of chunk glm-better

  • 観察データは常に正の値なのに予測が負に突入してない?
  • 縦軸は整数。しかものばらつきが横軸に応じて変化?
  • データに合わせた統計モデルを使うとマシ

ちょっとずつ線形モデルを発展させていく

線形モデル LM (単純な直線あてはめ)

    ↓ いろんな確率分布を扱いたい 👈 統計モデルの重要な部品

一般化線形モデル GLM

    ↓ 個体差などの変量効果を扱いたい

一般化線形混合モデル GLMM

    ↓ もっと自由なモデリングを!

階層ベイズモデル HBM

データ解析のための統計モデリング入門 久保拓弥 2012 より改変

確率分布

発生する事象(値)と頻度の関係。

手元のデータを数えて作るのが経験分布
e.g., サイコロを12回投げた結果、学生1000人の身長

plot of chunk distribution

一方、少数のパラメータと数式で作るのが理論分布
(こちらを単に「確率分布」と呼ぶことが多い印象)

確率変数$X$はパラメータ$\theta$の確率分布$f$に従う…?

$X \sim f(\theta)$

e.g.,
コインを3枚投げたうち表の出る枚数 $X$ は二項分布に従う
$X \sim \text{Binomial}(n = 3, p = 0.5)$

plot of chunk dbinom

\[\begin{split} \text{Prob}(X = k) &= \binom n k p^k (1 - p)^{n - k} \\ k &\in \{0, 1, 2, \ldots, n\} \end{split}\]

一緒に実験してみよう。

試行を繰り返して記録してみる

コインを3枚投げたうち表の出た枚数 $X$

試行1: → $X = 2$
試行2: 裏 裏 裏 → $X = 0$
試行3: 裏 裏 → $X = 1$ 続けて $2, 1, 3, 0, 2, \ldots$

plot of chunk rbinom

試行回数を増やすほど二項分布の形に近づく。
0と3はレア。1と2が3倍ほど出やすいらしい。

コイントスしなくても $X$ らしきものを生成できる

  • コインを3枚投げたうち表の出る枚数 $X$
  • $n = 3, p = 0.5$ の二項分布からサンプルする乱数 $X$
plot of chunk dbinom
$X \sim \text{Binomial}(n = 3, p = 0.5)$

   ↓ サンプル

{2, 0, 1, 2, 1, 3, 0, 2, …}

これらはとてもよく似ているので
「コインをn枚投げたうち表の出る枚数は二項分布に従う」
みたいな言い方をする。逆に言うと
「二項分布とはn回試行のうちの成功回数を確率変数とする分布」
のように理解できる。

統計モデリングの一環とも捉えられる

コイン3枚投げを繰り返して得たデータ {2, 0, 1, 2, 1, 3, 0, 2, …}

   ↓ たった2つのパラメータで記述。情報を圧縮。

$n = 3, p = 0.5$ の二項分布で説明・再現できるぞ


「データ分析のための数理モデル入門」江崎貴裕 2020 より改変

こういうふうに現象と対応した確率分布、ほかにもある?

有名な確率分布、それに「従う」もの

離散一様分布
コインの表裏、サイコロの出目1–6
幾何分布
成功率pの試行が初めて成功するまでの失敗回数
二項分布
成功率p、試行回数nのうちの成功回数
ポアソン分布
単位時間あたり平均$\lambda$回起こる事象の発生回数
ガンマ分布
ポアソン過程でk回起こるまでの待ち時間
(k = 1のとき指数分布と呼ばれる)
正規分布
確率変数の和、平均値

離散一様分布

同じ確率で起こるn通りの事象のうち実際に起こった事象X

e.g., コインの表裏、サイコロの出目1–6

plot of chunk dunif

🔰 一様分布になりそうな例を考えてみよう

幾何分布 $~\text{Geom}(p)$

成功率pの試行が初めて成功するまでの失敗回数X

e.g., コイントスで表が出るまでに何回裏が出るか

plot of chunk geometric

\[ \text{Prob}(X = k \mid p) = p (1 - p)^k \]

「初めて成功するまでの試行回数」とする定義もある。

🔰 幾何分布になりそうな例を考えてみよう

二項分布 $~\text{Binomial}(n,~p)$

確率$p$で当たるクジを$n$回引いたうち当たった回数X。平均は$np$。

plot of chunk dbinom-n

\[ \text{Prob}(X = k \mid n,~p) = \binom n k p^k (1 - p)^{n - k} \]

🔰 二項分布になりそうな例を考えてみよう

ポアソン分布 $~\text{Poisson}(\lambda)$

単位時間(空間)あたりに平均$\lambda$回発生する事象が実際に起きた回数X。

e.g., 1時間あたりのメッセージ受信件数、メッシュ区画内の生物個体数

plot of chunk dpoisson

\[ \text{Prob}(X = k \mid \lambda) = \frac {\lambda^k e^{-\lambda}} {k!} \]

二項分布の極限 $(\lambda = np;~n \to \infty;~p \to 0)$。
めったに起きないことを何回も試行するような感じ。

指数分布 $~\text{Exp}(\lambda)$

ポアソン過程の事象の発生間隔x。平均は $1 / \lambda$ 。

e.g., メッセージの受信間隔、道路沿いに落ちてる手袋の間隔

plot of chunk dexp

\[ \text{Prob}(x \mid \lambda) = \lambda e^{-\lambda x} \]

幾何分布の連続値版。

🔰 ポアソン分布・指数分布になりそうな例を考えてみよう

ガンマ分布 $~\text{Gamma}(k,~\lambda)$

ポアソン過程の事象k回発生までの待ち時間x

e.g., メッセージを2つ受信するまでの待ち時間

plot of chunk dgamma

\[ \text{Prob}(x \mid k,~\lambda) = \frac {\lambda^k x^{k - 1} e^{-\lambda x}} {\Gamma(k)} \]

指数分布をkのぶん右に膨らませた感じ。
shapeパラメータ $k = 1$ のとき指数分布と一致。

正規分布 $~\mathcal{N}(\mu,~\sigma)$

平均 $\mu$、標準偏差 $\sigma$ の美しい分布。よく登場する。
e.g., $\mu = 50, ~\sigma = 10$ (濃い灰色にデータの95%, 99%が含まれる):

plot of chunk gaussian

\[ \text{Prob}(x \mid \mu,~\sigma) = \frac 1 {\sqrt{2 \pi \sigma^2}} \exp \left(\frac {-(x - \mu)^2} {2\sigma^2} \right) \]

正規分布に近づくものがいろいろある

標本平均の反復(中心極限定理); e.g., 一様分布 [0, 100) から40サンプル

plot of chunk central-limit

大きい$n$の二項分布

plot of chunk binom-normal

正規分布に近づくものがいろいろある

大きい$\lambda$のポアソン分布

plot of chunk poisson-normal

平均値固定なら$k$が大きくなるほど左右対称に尖るガンマ分布

plot of chunk gamma-normal

有名な確率分布対応関係ふりかえり


離散一様分布
コインの表裏、サイコロの出目1–6
幾何分布
成功率pの試行が初めて成功するまでの失敗回数
二項分布
成功率p、試行回数nのうちの成功回数
ポアソン分布
単位時間あたり平均$\lambda$回起こる事象の発生回数
ガンマ分布
ポアソン過程でk回起こるまでの待ち時間
(k = 1のとき指数分布と呼ばれる)
正規分布
確率変数の和、平均値。使い勝手が良く、よく登場する。

現実には、確率分布に「従わない」ことが多い

植物100個体から8個ずつ種子を取って植えたら全体で半分ちょい発芽。
親1個体あたりの生存数はn=8の二項分布になるはずだけど、
極端な値(全部死亡、全部生存)が多かった。

plot of chunk overdispersion

「それはなぜ?」と考えて要因を探るのも統計モデリングの仕事。
「普通はこれに従うはず」を理解してこそできる思考。

疑似乱数生成器 Pseudo Random Number Generator

コンピューター上でランダムっぽい数値を出力する装置。
実際には決定論的に計算されているので、
シード(出発点)と呼び出し回数が同じなら出る数も同じになる。

set.seed(42)
runif(3L)
# 0.9148060 0.9370754 0.2861395
runif(3L)
# 0.8304476 0.6417455 0.5190959
set.seed(42)
runif(6L)
# 0.9148060 0.9370754 0.2861395 0.8304476 0.6417455 0.5190959

シードに適当な固定値を与えておくことで再現性を保てる。
ただし「このシードじゃないと良い結果が出ない」はダメ。

さまざまな「分布に従う」乱数を生成することもできる。

いろんな乱数を生成・可視化して感覚を掴もう

n = 100
x = sample.int(6, n, replace = TRUE)
x = runif(n, min = 0, max = 1)
x = rgeom(n, prob = 0.5)
x = rbinom(n, size = 3, prob = 0.5)
x = rpois(n, lambda = 10)
x = rnorm(n, mean = 50, sd = 10)
print(x)

p1 = ggplot(data.frame(x)) + aes(x)
p1 + geom_histogram() # for continuous values
p1 + geom_bar()       # for discrete values

🔰 小さい n から徐々に大きくして変化を確認しよう。

🔰 ほかのオプションもいろいろ変えて変化を確認しよう。

🔰 1%の当たりを狙って10連ガチャを回す人が100万人いたら、
全部はずれ、1つ当たり、2つ当たり… の人はどれくらいいるか?

データに分布をあてはめたい

ある植物を50個体調べて、それぞれの種子数Xを数えた。

plot of chunk poisson-seed

カウントデータだからポアソン分布っぽい。
ポアソン分布のパラメータ $\lambda$ はどう決める?

データに分布をあてはめたい

ある植物を50個体調べて、それぞれの種子数Xを数えた。

plot of chunk poisson-seed-lambda

カウントデータだからポアソン分布っぽい。
ポアソン分布のパラメータ $\lambda$ はどう決める?
(黒が観察データ。青がポアソン分布。よく重なるのは?)

ゆう (likelihood)

もっともらしさ。 モデルのあてはまりの良さの尺度のひとつ。

あるモデル$M$の下でそのデータ$D$が観察される確率
定義通り素直に書くと
$\text{Prob}(D \mid M)$

データ$D$を固定し、モデル$M$の関数とみなしたものが尤度関数:
$L(M \mid D)$

モデルの構造も固定してパラメータ$\theta$だけ動かす場合はこう書く:
$L(\theta \mid D)$ とか $L(\theta)$ とか

尤度を手計算できる例

コインを5枚投げた結果 $D$: 表 4, 裏 1

表が出る確率 $p = 0.5$ と仮定:

\[\begin{split} L(0.5 \mid D) &= \binom 5 1 \times \text{Prob}(表 \mid 0.5) ^ 4 \times \text{Prob}(裏 \mid 0.5) ^ 1 \\ &= 5 \times 0.5 ^ 4 \times 0.5 ^ 1 = 0.15625 \end{split}\]

表が出る確率 $p = 0.8$ と仮定:

\[\begin{split} L(0.8 \mid D) &= \binom 5 1 \times \text{Prob}(表 \mid 0.8) ^ 4 \times \text{Prob}(裏 \mid 0.8) ^ 1 \\ &= 5 \times 0.8 ^ 4 \times 0.2 ^ 1 = 0.4096 \end{split}\]

$L(0.8 \mid D) > L(0.5 \mid D)$

$p = 0.8$ のほうがより尤もらしい。

種子数ポアソン分布の例でも尤度を計算してみる

ある植物が作った種子を数える。$n = 50$個体ぶん。

\[\begin{split} L(\lambda \mid D) = \prod _i ^n \text{Prob}(X_i \mid \lambda) = \prod _i ^n \frac {\lambda ^ {X_i} e ^ {-\lambda}} {X_i !} \end{split}\]

plot of chunk poisson-seed-likelihood

この中では $\lambda = 3$ がいいけど、より尤もらしい値を求めたい。

最尤推定 Maximum Likelihood Estimation

扱いやすい 対数尤度 (log likelihood) にしてから計算する。
一階微分が0になる $\lambda$ を求めると…この場合は標本平均と一致。

\[\begin{split} \log L(\lambda \mid D) &= \sum _i ^n \left[ X_i \log (\lambda) - \lambda - \log (X_i !) \right] \\ \frac {\mathrm d \log L(\lambda \mid D)} {\mathrm d \lambda} &= \frac 1 \lambda \sum _i ^n X_i - n = 0 \\ \hat \lambda &= \frac 1 n \sum _i ^n X_i \end{split}\]

plot of chunk poisson-mle

最尤推定を使っても“真のλ”は得られない

今回のデータは真の生成ルール“$X \sim \text{Poisson}(\lambda = 3.0)$”で作った。
「50個体サンプル→最尤推定」を1,000回繰り返してみると:

plot of chunk poisson-mle-repl

サンプルの取れ方によってはかなりズレた推定をしてしまう。
(標本データへのあてはまりはかなり良く見えるのに!)

サンプルサイズを増やすほどマシにはなる

“$X \sim \text{Poisson}(\lambda = 3.0)$”からnサンプル→最尤推定を1,000回繰り返す:

plot of chunk poisson-mle-nsam

Q. じゃあどれくらいのサンプル数nを確保すればいいのか?
A. 推定したい統計量とか、許容できる誤差とかによる。

すべてのモデルは間違っている

確率分布がいい感じに最尤推定できたとしても、
それはあくまでモデル。仮定。近似。

All models are wrong, but some are useful. — George E. P. Box


「データ分析のための数理モデル入門」江崎貴裕 2020 より改変

統計モデリングの道具 — まとめ

  • 何はともあれ作図して俯瞰
  • 確率変数 $X$
  • 確率分布 $X \sim f(\theta)$
    • 少ないパラメータ $\theta$ でばらつきの様子を表現
    • この現象はこの分布を作りがち(〜に従う) という知見がある
  • 尤度
    • あるモデルでこのデータになる確率 $\text{Prob}(D \mid M)$
    • データ固定でモデル探索 → 尤度関数 $L(M \mid D),~L(\theta \mid D)$
    • 対数を取ったほうが扱いやすい → 対数尤度 $\log L(M \mid D)$
    • これを最大化するようなパラメータ $\hat \theta$ 探し = 最尤法

🔰 3日目の課題2

サイコロを10回振ったら6の目が3回出た。

  1. 6の目の出る確率が1/6だとした場合の尤度は?

  2. 6の目の出る確率が0.2だとした場合の尤度は?

  3. 横軸を6の目の出る確率、縦軸を対数尤度とするグラフを描こう。

  4. このサイコロで6の目が出る確率を最尤推定しよう。
    数学で解ければ優。Rで見つければ良。目分量・勘で可。

ヒント
確率pで当たるクジをn回引いてk回当たる確率、と同じ計算を使う。
数学の $\binom 5 2 = {}_5 \mathrm{C} _2 = 10$ はRでは choose(5, 2)

参考文献