2.1 例題:種子数の統計モデリング
setwd("~/Desktop/Statistics_Study/StatisticalModeling/kubobook_2012/02distribution")
load("data.RData")#種子数のデータ(1つの花が何個種子を持ってるか)
head(data)
## [1] 2 2 4 6 4 5
length(data)
## [1] 50
summary(data)#要約統計量
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 2.00 3.00 3.56 4.75 7.00
table(data)#度数分布
## data
## 0 1 2 3 4 5 6 7
## 1 3 11 12 10 5 4 4
par(family = "HiraKakuProN-W3") #日本語もOKに
hist(data, breaks=seq(-0.5, 9.5, 1), main="種子数のヒストグラム") #ヒストグラム(-0.5から9.5を1ずつ)
var(data) #分散
## [1] 2.986122
sd(data) #標準偏差
## [1] 1.72804
#ggplot2を使ってみた
library(ggplot2)
#dataframeへの変換が必要
dataframe<-data.frame(count = data)
ggplot(data = dataframe) +
aes(x = count) +
xlab("count_data") +
ylab("Frequency") +
geom_bar()
2.2 データと確率分布の対応関係を眺める
#平均3.56のポアソン分布
#種子数0~9
y <-0:9
prob <- dpois(y, lambda = 3.56)
cbind(y, prob)
## y prob
## [1,] 0 0.02843882
## [2,] 1 0.10124222
## [3,] 2 0.18021114
## [4,] 3 0.21385056
## [5,] 4 0.19032700
## [6,] 5 0.13551282
## [7,] 6 0.08040427
## [8,] 7 0.04089132
## [9,] 8 0.01819664
## [10,] 9 0.00719778
hist(data, col="white")
par(new=T)
#typeが点と線, ltyで破線に
plot(y, prob, type = "b", lty = 5, xlab="", ylab="", axes=F)
axis(4)
mtext("Frequency",side = 4,
#どこにラベルを置くか(1なら下,2なら左,3なら上,4なら右),
line = 2 #グラフの枠からの距離
)
2.3 ポアソン分布とは何か
#例えばラムダ=3.5
plot(0:20, dpois(0:20, 3.5), type="b", lty=5)
2.4 ポアソン分布のパラメータの最尤推定
#data=data, 平均mのポアソン分布(対数)の確率密度の合計
logL<- function(m) {
sum(dpois(data, m, log=TRUE))
}
lambda<- seq(2, 5, 0.1) #平均を2~5で0.1ずつ
#x=lambda sapplyでlambdaの各オブジェクトにlogLをあてはめ
plot(lambda, sapply(lambda, logL), type="l")
abline(v=3.56)
- ポイント
- 対数尤度が0に近いほどヒストグラムに近そう
- λ=3.56が最尤推定値
2.4.1 擬似乱数と最尤推定値のばらつき
#平均3.5のポアソン分布に従うような乱数を50個発生させる(結果は毎回変わる)
#変えたくない場合はset.seed(整数)
r_data<-rpois(50, 3.5)
mean(r_data)
## [1] 3.82
2.6 確率分布の選び方
- 分析の時はどの確率分布になりそうかを考える
- 離散型? 連続型?
- 範囲
- 標本分散と標本平均の関係
2.7 この章のまとめと参考文献
- どういう理由でこの統計モデルを使いましたと言えるようになろう