ガンマ分布のscaleとshape(rate)を直感的にわかりやすく理解する【R】
統計モデリングを学ぶ上での名著、
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
(緑本)の6章8節にこのような記述がある。
この例題は、架空植物50個の葉の重量と花の重量の関係を調べます。個体
の葉重量を
、花重量を
とするのですが、
ここでは、ある個体の花の重量がガンマ分布に従っていることにします。応答変数は連続地ですが、重量なので正の値しかとりません。したがって、そのばらつきは正規分布ではなくガンマ分布で説明したほうがよさそうです。
なるほど、応答変数が連続値で正の値しかとらない場合はガンマ分布のほうがふさわしいのか。そのガンマ分布とやらを一度Rで見てみようということで、正規分布のグラフを書くノリでcurve()+dgammaを渡してみたところ、
curve(dgamma,0,10)
## Error in dgamma(x): 引数 "shape" がありませんし、省略時既定値もありません
当然ながら怒られました。ガンマ分布をR上で使う場合、shapeとscaleというパラメータを設定してあげる必要があります。じゃあそのshapeとscaleって何？というお話。（scaleの代わりにrateを渡すこともできます。scale=1/rateで、私はscaleの方が理解しやすいのでそっちで説明していきます。）
文字通り、「形」と「縮尺」
shapeとは形状母数でscaleとは尺度母数！
いってしまえばこれだけなのですが、細かく説明していきます。ただしガンマ分布の密度関数など、数式の話は一切なしです。ひたすら直感的に。
shape=どれだけ正規分布に近づけるか
ガンマ分布は先ほども述べたように、「正の値しかとらない連続値の分布」です。左端が決まっている分、なんとなく左右非対称(左に偏った形状)になっていることは想像できます。
その「左右非対称さ」を決めるのがshapeで、大きければ大きいほど左右対称(正規分布)に近づいていきます。小さいと左に偏っていきます(最小は0、直線になる)。
scaleは5に固定してshapeを1,2,3,4と上げていきましょう。
curve(dgamma(x , shape = 1 , scale = 5) , 1 , 30 , col = "black") curve(dgamma(x , shape = 2 , scale = 5) , 2 , 30 , add = T , col = "red") curve(dgamma(x , shape = 3 , scale = 5) , 1 , 30 , add = T , col = "blue") curve(dgamma(x , shape = 4 , scale = 5) , 1 , 30 , add = T , col = "green")
徐々に左右対称に近づいていっているのがわかります。軸を一緒にして重ねると何となく形の特徴が捉えにくいので別バージョンも。今度はshapeの値は左上から1,3,10,50にしています。
par(mfrow=c(2,2)) curve(dgamma(x , shape = 1 , scale = 5) , 1 , 30 , col = "black") curve(dgamma(x , shape = 3 , scale = 5) , 1 , 100 , col = "red") curve(dgamma(x , shape = 10 , scale = 5) , 1 , 150 , col = "blue") curve(dgamma(x , shape = 50 , scale = 5) , 100 , 400 , col = "green")
だいぶ変化がわかりやすくなりました。shape=50のグラフはもうほぼ平均250,分散1250の正規分布ですね(伏線)。
scale=目盛をどれだけ大きくするか
これはもっと簡単です。さっそく図を並べてみましょう。
shapeは5に固定、scaleは1,10,100,1000です。軸に注目。
par(mfrow=c(2,2)) curve(dgamma(x , shape = 5 , scale = 1) , 1 , 15 , col = "black") curve(dgamma(x , shape = 5 , scale = 10) , 1 , 150 , col = "red") curve(dgamma(x , shape = 5 , scale = 100) , 1 , 1500 , col = "blue") curve(dgamma(x , shape = 5 , scale = 1000) , 1 , 15000 , col = "green")
scaleがデカいほど目盛がデカい！単純！
じゃあ値に対してどのぐらい目盛がデカくなっているのかが気になるところです。
そう、今まで秘密にしてきましたがこのような関係がガンマ分布にはあるのです。なぜかというと、そういう風に定義されているからです(暴論)。
なのでどれくらい目盛がデカくなるかの答えは、「平均がshape×scaleになるぐらいの大きさになる」となります。
最初に私がrateよりscale派だといったのもそのためで、だいたい分布決めるときって平均考えながら決めるので、簡単にshapeと掛け合わせるだけで平均が求められるscaleのほうが何かと便利なのです。逆数ってイメージしずらくない？？
最後に実例を
今手元に、正の値をとる応答変数の、平均が100で若干左に偏った分布があり、ガンマ分布と比較してみようと思った時、どのようなプロットをしますか。
## ちょい左に偏ってるし`shape`は4ぐらいかな、なら平均は100だから`scale`は25になるな curve(dgamma(x, shape = 4 ,scale = 25) , 1 , 300)
