ガンマ分布のはなし - 統計学といくつかのよしなしごと

いまいちイメージがつかみにくい*1

　「ガンマ分布！」と言われてもイマイチなんのこっちゃわかりにくい気がする。二項分布やポアソン分布のように直感的なイメージを持ちにくく、教科書でも「ガンマ分布は指数分布の一般化です」などとどことなく味気ないからだと思う。なので今回は、ガンマ分布の成り立ちとこれに従う現実の現象をみることで、なんとなくイメージを持ちやすくしてもらうことを目指して記事を書いてみる。

ガンマ分布

　先に書いた通りガンマ分布は指数分布を一般化したもので、以下の確率密度関数で定義される。

${\displaystyle f(x) = \frac{λ^α}{Γ(α)} x^{α-1} e^{-λx} ~~ (x\geq0), ~~~~~ 0~~(x<0) }$

このガンマ分布に対して与えられる意味合いは、「一定期間に1回起きると期待されるランダムな事象が複数回起きるまでの時間の分布」というものだが、なぜこの式がその意味を持つのかを理解するには、幾何分布や指数分布との関連を見ていくのが良い。図示するとこんな感じ。赤枠の確率分布について、以下順に追ってみる。
f:id:yanbow221:20170902095612p:plain

幾何分布について

　幾何分布は成功（確率p）か失敗（確率q）のいずれかをとる試行（ベルヌーイ試行という）において、初めて成功するまでの試行回数の分布と定義される。その確率は、

${\displaystyle f(x) = pq^{x-1}, ~~x=1,2,3,... }$

初めて成功するまでの試行回数がx回なのだから(x-1)回の失敗と最後1回の成功ということでこれはわかりやすいと思う。この幾何分布は『離散的な待ち時間分布』とも呼ばれる。1回の成功を得る確率がpなのだから、別の見方をすれば1回の成功を得るためには(平均的に)1/pだけの時間を要するということからこの名前がついている。

指数分布について

　先ほどの幾何分布は試行回数と説明したように離散的な分布になっている。これを連続値について扱ったものが指数分布である。確率密度関数は、

${\displaystyle f(x) = λe^{-λx} ~~(x\geq0), ~~~~~ 0~~(x<0) }$

で定義される。幾何分布が離散的な待ち時間分布と呼ばれたのに対し、この指数分布は『連続的な待ち時間分布』と呼ばれる。幾何分布と指数分布は一見したときの式の形が異なっているため、どちらも待ち時間分布という同じ性質を持っていることを理解するには多少の証明が必要になる。以下がその解説。

<指数分布が連続的な待ち時間分布と呼ばれることのワケ>

実数xに対して、y=[x]をxを超えない最大の整数と定義する。例えばx=3.2ならy=3、x=5.84ならy=5、といった感じ。指数分布に従う確率変数Xに対して確率変数Yはどのような分布に従うかを考える。nを正の整数として、

${\displaystyle \begin{eqnarray} P(Y = n) &=& P([X ] = n) \\ &=& P(n \leq X < n+1) \end{eqnarray} }$

これより、

${\displaystyle \begin{eqnarray} P(Y \geq n) &=& P(Y = n) + P(Y = n+1) + ~~... \\ &=& P(n \leq X < n+1) + P(n+1 \leq X < n+2) + ~~... \\ &=& P(X \geq n) \\ &=& \int_n^∞ λe^{-λx} dx \\ &=& e^{-λn} \end{eqnarray} }$

よって、
${\displaystyle \begin{eqnarray} P(Y=n) &=& P(Y \geq n) - P(Y \geq n+1) \\ &=& e^{-λn} - e^{-λ(n+1)} \\ &=& e^{-λn} (1-e^{-λ}) \\ &=& (e^{-λ})^n (1-e^{-λ}) \end{eqnarray} }$

${\displaystyle q = e^{-λ} }$ とすれば、幾何分布の式と同じ形をしていることがわかる。この指数分布は、「期間1/λあたりに1回起きると期待される事象が1回起きるまでの時間の分布」という意味合いを持つことになる*2。

ガンマ分布について

　ガンマ分布は指数分布の一般化である。改めてガンマ分布の確率密度関数は以下。

${\displaystyle f(x) = \frac{λ^α}{Γ(α)} x^{α-1} e^{-λx} ~~ (x\geq0), ~~~~~ 0~~(x<0) }$

指数分布と見比べてみるとよく似た形をしていることがわかる。実際α=1としてやると指数分布に一致する。
※補足すると、分母のΓ(α)は階乗n!を一般化したもので、

${\displaystyle Γ(1)=1 }$

となる。このガンマ関数が式内に入っている理由は、確率密度関数を

${\displaystyle f(x) = \int_0^∞ \left \{ \frac{λ^α}{Γ(α)} x^{α-1} e^{-λx}\right \} dx = 1 }$

と1にするため(確率なので定義域全体で積分すると1になる必要がある)に必要な係数だと理解しておいていい。なお、詳しい導出方法は以下が参考になる。

ガンマ関数の導出

　さて、指数分布を一般化したものであるガンマ分布が持つ意味合いとしては、「期間1/λあたりに1回起こると期待されるランダムな事象がα回起こるまでの時間の分布」となる。ガンマ分布に従う現象としては、製品部品の寿命やウイルスの潜伏期間のほか、人の体重なんかも従うらしい*3。

　ここまで、ガンマ分布とその他の確率分布の関係性やその意味合いについて何となく見てきた。

ガンマ分布の例

　ここからは実際にデータを使ってガンマ分布がどんな形かを確認してみる。ただ少し探したところ適当なデータセットが見つからなかったので、今回はガンマ分布の意味合いから考えて闇のゲーム*4を題材にしてこれを考えることにする。

<闇のゲーム「モンスター・ワールド」とは>

　今回取り扱うのは、僕世代の男性なら知らない人はいないであろう『遊戯王』で、千年パズルを持つ主人公遊戯と千年リングを持つライバル獏良が戦ったボードゲームである。詳しくは単行本を読んでいただくとして、ここで知っておくべきルールは、

・獏良は1体のボス『ゾーク』を操り、遊戯は5体のキャラクターを操る
・遊戯は盤上でキャラクターを進めていってゾークを倒せば勝ち。逆にキャラクター全員が倒れると負け
・キャラクターとゾークは互いに10面ダイス2個を振って出た目の数で勝負をする。赤いダイスは十の位、白いダイスは一の位の数字を表す(つまり目は00~99の値をとる)。
・ダイスの目の数が「小さい」ほうが良い目で、良い目を出したほうが相手へ攻撃することができる(同じ目の時はやり直し)
・相手より良い目を出しかつその目が「00~09」だった場合、クリティカルヒットになる
・1度クリティカルヒットを出したキャラクターは勝ち抜けとなり、次のキャラクターとゾークのバトルに進む
・遊戯はクリティカルヒットを全キャラクターで出せば(つまり累計5回のクリティカルヒット)、ゾークを倒すことができ勝ちとなる
・獏良は『洗脳ダイス』により常に「10」の目を出してくる
・獏良のほうが良い目であっても、遊戯側のキャラクター『杏子(アンズ)』が即座に回復魔法をかけるため実質遊戯のキャラクターを倒すことはない
・単位時間(10秒)あたりのダイスロールは平均的に0.7回

である。図化するとこんな感じ。なお、ガンマ分布のシミュレーションのため原作のルールを多少(というか大幅に)変更している。
f:id:yanbow221:20170905235554p:plain

　このルールのもと、遊戯は獏良に勝つまでに大体どの程度の時間がかかるかを見てみる。まずキャラクターがゾークにクリティカルヒットを出すには「00, 01, 02, 03, ..., 09」のいずれかの目を出さなければならないので、その確率は1/10である。単位時間10秒あたりダイスロールは0.7回のペースで行われるので、単位時間当たりのクリティカルヒットの可能性は0.7×(1/10)=0.07、すなわちλ=0.07である。これは言い換えればキャラクター1人がクリティカルヒットを出すまでにかかる平均的な時間は(1/0.07)×10≒143秒ということである。これが各キャラクターにおいて出れば遊戯の勝ちだから、その時間の分布はλ=0.07、α=5のガンマ分布に従うと予想される。その確率密度関数は以下のようにして描かれる。

#闇のゲームのシミュレーション
set.seed(221)

#λ=0.035、α=5のガンマ分布
x <- seq(0, 600, 1)
λ <- 0.07
α <- 5

curve(dgamma(x, shape = α, scale = 1/λ),
      xlim = c(0, 600), ylim = c(0, 0.02), xlab = "時間(×10秒)", ylab = "確率密度")
par(new = T)
abline(v = 71.4, col = "red") #平均値

f:id:yanbow221:20170905232317j:plain

これに対して10,000回シミュレーションを回してみてその分布を比較してみるとこんな感じ。

#シミュレーション
hako <- matrix(nrow = 10000, ncol = 1) #結果の格納箱
dice <- seq(0, 99, 1)

for(i in 1:10000){
  count <- 0 #クリティカルヒット数カウント用
  total <- 0 #全ダイスロール回数カウント用
  while(count <= 5){
    x <- sample(dice, 1)
    if(x < 10){
      count <- count + 1
      total <- total + 1
    }else{
      total <- total + 1
    }
  }
  hako[i, 1] <- total
}

hako2 <- (hako - 5) / 0.7 #実質離散分布で、5単位未満で遊戯が勝つ可能性はないため調整

dens <- density(hako2) #確率密度の推定
par(new = T)
lines(dens, col = "blue", xlim = c(0, 600)) #シミュレーションに基づく分布

f:id:yanbow221:20170905233156j:plain

今回はシミュレーションの方法的に実質離散型の分布になっているため(だと思う)多少ずれているが、概ねフィットしているといえそう。遊戯は獏良から勝利を得るのに大体12分ほどかかっている。ちなみに得られたシミュレーションデータについてガンマ分布を仮定したパラメータ推定を行うとα'=5.6、λ'=0.07となる。

> est <- fitdistr(hako2, "gamma")
> est
      shape         rate    
  5.615397633   0.071035304 
 (0.077115818) (0.001020378)

終わりに

　今回は式が複雑でややイメージを掴みにくいガンマ分布についてみてきた。このガンマ分布はベイズ統計でも活躍する分布なので、こちらの勉強がある程度形になったらまた取り上げたい。
　ところで、冒頭で分布の関連を図示した中で指数分布とポアソン分布の間に双方向の矢印を引いた。実はこの2つの分布は同じ事象を別の観点から見たものになる。これについては以下の資料がとても分かりやすい。

指数分布とポアソン分布のいけない関係 from tera monagi

www.slideshare.net

　ガンマ分布について見たので次はベータ分布についてなんか書いてみよう。

*1:個人の感想です。

*2:確率変数Xの期待値E(X)を求めることで期間1/λというのが導かれる

*3:体重については待ち時間分布としての性質というより、(1)一般に身長が正規分布に従うこと、(2)ガンマ分布の三乗根は正規分布になること、(3)人の縦横高さの積で体重を近似できそうなこと、から理解できるようなのだけど、この(2)について詳しく確認できていない。詳しい方いらっしゃれば、教えていただけると嬉しいです。m(_ _)m

*4:小学生の頃はこれとMTGというカードゲームになけなしのお小遣いをつぎ込んだ