みどり本攻略 Day1
Rの復習と統計の勉強を兼ねて.みどり本とは「データ解析のための統計モデリング入門」という素晴らしい本です.学部時代にざっと読んだのですが,結局字面を追うだけでよく理解していなかったのでこの際再度しっかり学ぶことにしました.1日1章くらいのペースで読んでいければいいなと思います.基本的には自分のノートとして記述しますが,間違いなどあれば指摘していただけると幸いです.
この本の読者として想定しているのは...「数理モデルで現象を表現・説明する」基礎訓練を受けていない人たちです.
（まえがき より）
まさしく私のことです.
人は自然現象を観測しますが,その観測の数や程度は限られています.だから観測データをもとにその裏に隠された自然現象の法則（まさしく神さまのデザイン）を明らかにするために統計モデルを作って普遍化することが必要なのかなと思います.医学の世界では結構適当な統計もどきで有意差が出て大喜びということがありうるので,自らを律する意味でもしっかり学びたいと思います.
***
さて,初日の今日は２章を読みます.もちろんそのままコピーするのは良くないので,自分なりに学んだポイントを挙げていきます.
みどり本ではある植物の種子の数をもとに統計モデルを作るという題材で,解説していますが,最初のサンプルデータを変えて挑戦してみます.ある植物50個体の種子の数の架空データを使用しますが,あえて100個体分のデータを作って同じようにやってみます.
ネタバレ感がありますが,最初にデータを生成します.
data <- rpois(100, lambda=4.5)
これで data にポアソン分布に従う平均4.5の擬似乱数が100個格納されました.
このdataについてRの関数を使用して,さまざまな統計量を調べていきます.
summary(data)
Rはsummary()で色々出してくれるのがありがたいですね.
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 3.00 5.00 4.77 6.00 14.00
次に平均値を出してみます.
mean(data)
[1] 4.77
今回のデータは平均値が4.77です.
標本統計量も出してみます.
標本分散 sample variance
var(data)
[1] 5.512222
標本分散は,標本の相加平均をそれぞれの標本から引いて二乗した値の相加平均.
その平方根が標本標準偏差 sample standard deviation
sd(data)
sqrt(var(data)) #sqrt()は平方根を計算する
[1] 2.347812
さてここまででわかった今回のデータn=100について
平均値（相加平均）：4.77
中央値：5.00
最大値：14.00
最小値：0.00
25パーセンタイル値：3.00
75パーセンタイル値：6.00
標本分散：5.512222
標本標準偏差：2.347812
以上のことがわかりました.
次にこのデータを統計モデルに当てはめていきます.
さてみどり本の記述では
統計学における推定.自然を確率分布であると見立てて,得られたデータをそこから発生した乱数のセットだと考える.限定された個数のデータから元の確率分布に近い確率分布を探すのが「推定」
とあります.今回のデータを「自然」としてこれを確率分布であると見立てます.
さて確率分布 probability distributionとは確率変数 random variableの値と出現する確率を対応させたもの,です.確率分布にはいろいろなものがあり,それぞれにパラメータが設定されています.ここではポアソン分布を使用します.確率変数はここで言うところのある個体の種子数といったそれぞれにばらつく値です.
ポアソン分布 Poisson distribution
所与の時間間隔で発生する離散的な事象を数える特定の確率変数 X を持つ離散確率分布のことである。(Wikipediaより)
意味がわからないので、具体的な例を挙げると
ということがあるようです.
ポアソン分布を用いる状況は,
(1) データが正の整数と0である. 離散的.
(2) 下限はわかっているが上限はわからない.
(3) データの平均と分散が大体等しい.
今回のでデータは平均が4.77なので平均4.77のポアソン分布をグラフとして図示してみようと思います.この流れ自体はみどり本に詳しい記載があります.
y <- 0:14
prob <- dpois(y, lambda = 4.77)
これでprobには種子の数が0個から14個それぞれの場合の平均4.77のポアソン分布に従う確率が格納されます.
prob
[1] 0.0084803802 0.0404514134 0.0964766209 0.1533978272 0.1829269089
[6] 0.1745122711 0.1387372555 0.0945395298 0.0563691947 0.0298756732
[11] 0.0142506961 0.0061796200 0.0024563990 0.0009013095 0.0003070890
こんな感じです.さて図示してみます.
plot(y, prob, type="b", lty=2) #"b"は丸と折れ線,2は折れ線を破線にする
できた図はこれです.
これを見ると一番出現確率が高そうなのは4でしょうか.
では確率分布のパラメータを観測データに基づいて推定する方法として最尤推定を考えます.尤度の概念はすごくわかりやすい解説を発見したので下記に記します.
こちらを発見.
尤度関数の基本概念は、「サンプリングしてデータが観測された後、そのデータは元々どういうパラメーターを持つ確率分布から生まれたものだったか？」と言う問いに答えるためのものです。
その後も非常にわかりやすい解説が並びます.勉強になりました.
とりあえず尤度自体は扱いづらい数（すごく小さくなって小数点以下が長くなる）なので,対数変換して対数尤度関数 log likelihood function を使います.
尤度は実際のところ,それぞれの確率変数が出現する確率の積です.
さて,みどり本に従ってRを動かします.
logL <- function(m) sum(dpois(data, m, log =TRUE))
まず関数を定義します.Rで関数を定義するときはfunction(変数)の後にその内容を記述します.dpois(data,m,log=TRUE)では任意のパラメータmのポアソン分布において,dataを確率変数としたときの確率密度が出力されます.この際log=TRUEとしているので確率は対数で表されます.尤度は確率の積ですから,logにすると足し算となってsum()を使用することとなります.
lambda <- seq(2, 7, 0.1)
次にパラメータとなるlambdaを設定します.今回は大体4~5くらいで最大になりそうなので,2から7の間で0.1刻みの数を代入します.
plot(lambda, sapply(lambda, logL), type="l"
sapply()はある数列に対して,特定の関数を適用する関数です.ここではlambdaに格納されている数をlogLの変数に代入した結果をy軸に,lambdaをx軸にプロットした図が出力されます.
さて,このときの最尤推定値はこのグラフの傾きが0となるところですから,対数尤度関数を偏微分してイコールゼロとなるときのlambdaを求めるわけですが,この辺はみどり本通りにデータの標本平均となり,4.77となります.
大体こんな感じでしょうか.やっと一歩踏み出した感じですが,引き続き頑張ります.
