In [1]:
date()
In [2]:
sapply(c("pipeR", "ggplot2", "dplyr", "tidyr", "readr"), require, character.only = TRUE)
In [3]:
options(repr.plot.width = 4, repr.plot.height = 4)
2章¶
種子数の統計モデリング¶
例題
- 調査対象種の植物の複数個体から種子を採取
- 個体$i$から$y_i$個の種子を採取
- $y_i$はどのような確率分布に従うか
例題データ読み込み（http://hosho.ees.hokudai.ac.jp/%7Ekubo/stat/iwanamibook/fig/distribution/data.RData）
In [4]:
dir(".", full.names = TRUE) %>>% cat(sep = "\n")
In [5]:
dir("./data", full.names = TRUE, recursive = TRUE) %>>% cat(sep = "\n")
In [6]:
ls.str()
In [7]:
load("data/chap02/data.RData")
In [8]:
ls.str()
データはすべて非負の整数
In [9]:
data
個数は50（50個体分のデータ）
In [10]:
length(data)
標本平均は3.56
In [11]:
summary(data)
度数分布
In [12]:
table(data)
ヒストグラムを描く
- ひと山の分布
In [13]:
data_frame(data = data) %>>%
ggplot(aes(x = data)) +
geom_histogram(binwidth = 1, colour = "black", fill = gray(0.6)) +
scale_y_continuous(breaks = seq(0, 12, 2)) +
scale_x_continuous(breaks = seq(0, 8, 2), limits = c(-0.5, 8.1)) +
labs(title = "Histogram of data") +
ylab("Frequency")
標本分散
In [14]:
var(data)
標本標準偏差
In [15]:
sd(data)
In [16]:
var(data) %>>% sqrt
2.2 データと確率分布の対応関係をながめる¶
- 確率分布: 確率変数の値とそれが出現する確率を対応させたもの
- 確率変数: ある植物個体 $i$ の種子数 $y_i$
- 個体$i$の種子数が$y_i$になる確率はポアソン分布に従う，と考える
- ポアソン分布]
- ポアソン分布 - Wikipedia
- パラメータ: 分布の平均
平均 3.56 のポアソン分布は，
In [17]:
y <- 0:9
prob <- dpois(y, lambda = 3.56)
- 種子数が0になる確率は，0.028，1になる確率は0.10
- 一番確立が高いのは種子数が3になる場合で，0.21
In [18]:
cbind(y, prob)
In [19]:
data_frame(y, prob) %>>%
ggplot(aes(y, prob)) +
geom_line(linetype = 2) +
geom_point() +
scale_x_continuous(breaks = seq(0, 8, 2))
観測データのヒストグラムと確率分布を重ねる（確率分布 $\times$ 個体数）
In [20]:
ggplot() +
geom_histogram(data = data_frame(x = data), aes(x = x),
colour = "black", fill = gray(0.6), binwidth = 1) +
geom_line(data = data_frame(y, prob = prob * 50), aes(x = y, y = prob), linetype = 2) +
geom_point(data = data_frame(y, prob = prob * 50), aes(x = y, y = prob),
shape = 21, fill = "white") +
labs(title = "Histogram of data") +
ylab("Frequency") +
xlab("data") +
scale_x_continuous(breaks = seq(0, 8, 2)) +
scale_y_continuous(breaks = seq(0, 12, 2))
2.3 ポアソン分布とは何か?¶
定義は，
$$p(y | \lambda) = \frac{\lambda^y \exp(-\lambda)}{y!} $$性質は，
- $\sum\limits^{\infty}_{y = 0}p(y | \lambda) = 1$
- 平均は $\lambda$
- 分散と平均は等しい
種子数のデータを表現するのにポアソン分布を選んだのは，
- 値が非負の整数
- $y_i$ に下限（0）はあるが上限はわからない
- 平均（3.56）と分散（2.99）がだいたい等しい
パラメータ $\lambda$ を変化させると，
In [21]:
data_frame(x = c(0:20), y1 = dpois(x, lambda = 3.5), y2 = dpois(x, lambda = 7.7),
y3 = dpois(x, lambda = 15.1)) %>>%
gather(lambda, val, -x) %>>%
ggplot(aes(x = x, y = val, shape = lambda)) +
geom_line() +
geom_point(fill = "white") +
ylab("prob") +
xlab("y") +
scale_shape_manual(values = c(21, 23, 24), labels = c("3.5", "7.7", "15.1")) +
theme(legend.position=c(.9, .85))
2.4 ポアソン分布のパラメータの最尤推定¶
- 最尤推定
- 尤度（当てはまりの良さ）を最大にするパラメータ（ここでは$\lambda$）の値を探す
- 尤度（$L(\lambda)$）はある$\lambda$ における $p(y_i | \lambda)$ の積 （同時確率）
- 実際には対数尤度関数を使う
$y_i$，$y_2$，$y_3$についてなら，尤度は，
In [22]:
data[1:3]
In [23]:
dpois(data[1:3], lambda = 3.56)
In [24]:
dpois(data[1:3], lambda = 3.56) %>>% prod()
平均 $\lambda$ を変化させたときにポアソン分布の形状と対数尤度がどのように変化するか
対数尤度にして計算
In [25]:
logL <- function(m){
dpois(data, m, log = TRUE) %>>% sum
}
In [26]:
gpPois <- function(lambda){
gp <- ggplot() +
geom_histogram(data = data_frame(x = data), aes(x = x),
colour = "black", fill = gray(0.6), binwidth = 1) +
geom_line(data = data_frame(y, prob = dpois(y, lambda = lambda) * 50), aes(x = y, y = prob),
linetype = 2) +
geom_point(data = data_frame(y, prob = dpois(y, lambda = lambda) * 50), aes(x = y, y = prob),
shape = 21, fill = "white") +
scale_x_continuous(breaks = seq(0, 8, 2)) +
scale_y_continuous(breaks = seq(0, 15, 5), limits = c(0, 15)) +
annotate(x = 6, y = 13, hjust = 0, vjust = 0, geom = "text", size = 3,
label = sprintf("lambda=%.1f\nlog L=%.1f", lambda, logL(lambda))) +
theme(axis.title = element_blank())
return(gp)
}
In [27]:
gps <- seq(2.0, 5.2, 0.4) %>>% lapply(gpPois)
In [28]:
library(gridExtra)
In [29]:
options(repr.plot.width = 8, repr.plot.height = 8)
$\lambda = 3.6$ のとき $\log L = -97.3$ で最大
In [30]:
do.call(grid.arrange, c(gps, ncol = 3))
対数尤度と$\lambda$の関係
- 3.5ぐらいで最大
In [31]:
options(repr.plot.width = 4, repr.plot.height = 4)
In [32]:
lambda = seq(2, 5, 0.1)
data_frame(lambda, logL = sapply(lambda, logL)) %>>%
ggplot(aes(x = lambda, y = logL)) +
geom_line() +
geom_vline(xintercept = 3.56, linetype = 2) +
scale_x_continuous(breaks = seq(2, 5, 0.5))
対数尤度が最大となる$\lambda$ を $\hat\lambda$ とする． 対数尤度関数が最大値で傾きゼロとなるとなる$\lambda$を探せばよい
$$ \begin{align} \frac{\partial \log L(\lambda)}{\partial \lambda} &= \sum\limits_{i}\left\{ \frac{y_i}{\lambda} - 1 \right\} \\ &= \frac{1}{\lambda}\sum\limits_{i}y_i - 50 \end{align} $$これが 0になるときの $\lambda$ が $\hat\lambda$ なので，解くと， $$ \hat\lambda = \frac{1}{50}\sum\limits_{i}y_i $$ これは標本平均のことなので $\hat\lambda = 3.56$ となる
- $\hat\lambda$: 最尤推定量
- $\hat\lambda = 3.56$: 最尤推定値
一般化すると，パラメータ$\theta$の確率分布から観測データ$y_i$が発生した時の確率を$p(y_i | \theta)$とすると尤度と対数尤度は，
$$ L(\theta | \boldsymbol{Y}) = \prod\limits_i p(y_i | \theta) $$$$ \log L(\theta | \boldsymbol{Y}) = \sum\limits_i \log p(y_i | \theta) $$となり，最尤推定で尤度（対数尤度）最大の$\hat{\theta}$を探す
In [33]:
binom.test(4, 10)
この区間の中で，尤度が最大となる$\theta$を求める（最尤推定）．
対数尤度は次のようになる．
$$\log L(\theta) = \log({}_{10}\mathrm{C}_4) + 4\log\theta + 6\log(1 - \theta)$$これをプロットすると（定数項は除く）
In [34]:
logL <- function(t){
4 * log(t) + 6 * log(1 - t)
}
ggplot(data_frame(x = c(0, 1)), aes(x)) +
stat_function(fun = logL) +
xlab(expression(theta)) +
ylab(expression("L("~theta~")")) +
geom_vline(xintercept = c(0.12, 0.74), linetype = 2)
$L(\theta)$ が最大になる$\theta$を求める．
$$ \frac{d}{d\theta} \log L(\theta) = \frac{4}{\theta} - \frac{6}{1-\theta} $$が0になる点なので，$\theta = 0.4$
（10枚中4枚表が出たので，表が出る確率は0.4）
In [35]:
data_frame(y = sapply(c(1:3000), function(x){
rpois(50, lambda = 3.5) %>>% mean
})) %>>% ggplot(aes(x = y)) +
geom_histogram(colour = "black", fill = gray(0.6), binwidth = 0.1) +
scale_x_continuous(limits = c(2.5, 4.5), breaks = seq(2.5, 4.5, 0.5)) +
xlab(expression("Estimated "~lambda~"for each trial"))
調査個体数が多ければ標準誤差は小さくなる
- 500個体
In [36]:
data_frame(y = sapply(c(1:3000), function(x){
rpois(500, lambda = 3.5) %>>% mean
})) %>>% ggplot(aes(x = y)) +
geom_histogram(colour = "black", fill = gray(0.6), binwidth = 0.1) +
scale_x_continuous(limits = c(2.5, 4.5), breaks = seq(2.5, 4.5, 0.5)) +
xlab(expression("Estimated "~lambda~"for each trial"))
2.5 統計モデルの要点: 乱数発生・推定・予測¶
- 観測データの裏には「真の統計モデル」がある
- 「真の統計モデル」はあるパラメータに従う確率分布
- その確率分布から生成された乱数を観測個数サンプリングしたものが観測データ
- パラメータを推定する or モデルをデータに当てはめる
- 推定した統計モデルを使って未知データを予測する
In [37]:
ggplot(data_frame(x = c(0:40)), aes(x)) +
stat_function(geom="point", n=41, fun = dbinom, args = list(size = 10, prob = 0.4), colour = "blue") +
stat_function(geom="point", n=41, fun = dbinom, args = list(size = 20, prob = 0.7), colour = "green") +
stat_function(geom="point", n=41, fun = dbinom, args = list(size = 40, prob = 0.5), colour = "red")
- 正規分布
- ガンマ分布
In [38]:
ggplot(data_frame(x = c(0, 20)), aes(x)) +
stat_function(fun = dgamma, args = list(shape = 1, rate = 0.5), colour = "red") +
stat_function(fun = dgamma, args = list(shape = 2, rate = 0.5), colour = "orange") +
stat_function(fun = dgamma, args = list(shape = 1, rate = 1), colour = "blue") +
stat_function(fun = dgamma, args = list(shape = 2, rate = 1), colour = "green")
- 一様分布
2.7 まとめ¶
- 確率分布で考える
- データへのモデルの当てはまりは尤度で評価できる
- 最尤推定: 尤度最大となるパラメータを推定すること
- 推定結果を用いた未知データへのあてはめ: 予測
- 確率分布を混合することで複雑なばらつきを表現できる
In [39]:
devtools::session_info()