データ解析のための統計モデリング入門 Chapter 2

参考文献

• http://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiBook.html

データ読み込み

load("~/statistics/modeling.for.data.analysis/data.RData")

データをみる

plant.seed <- data
length(plant.seed)

[1] 50

summary(plant.seed)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00    2.00    3.00    3.56    4.75    7.00 

hist(plant.seed, breaks = seq(-0.5, 9.5, 1))

var(plant.seed)

[1] 2.986

sd(plant.seed)

[1] 1.728

sqrt(var(plant.seed))

[1] 1.728

平均と分散が3前後で比較的近いので、ポアソン分布で近似できそう。

ポアソン分布をプロット

y <- 0:9
prob <- dpois(y, lambda = 3.56)
plot(y, prob, type = "b", lty = 2)

データとポアソン分布の対応をみる

hist(plant.seed, breaks = seq(-0.5, 9.5, 1))
points(y, prob * 50)
lines(y, prob * 50, lty = 2)

\( \lambda \)の値ごとのlog likelihoodをみる

logL <- function(m) sum(dpois(data, m, log = TRUE))

plot.poisson <- function(lambda) {
    y <- 0:9
    prob <- dpois(y, lambda = lambda)

    hist(plant.seed, breaks = seq(-0.5, 9.5, 1), ylim = c(0, 15),
         main = "", xlab = "", ylab = "")
    points(y, prob * 50)
    lines(y,  prob * 50, lty = 2)

    title(sprintf("lambda= %.1f\n logL= %.1f", lambda, logL(lambda)))
}

layout(matrix(1:9, byrow = T, ncol = 3))
junk <- sapply(seq(2, 5.2, 0.4), plot.poisson)

Likelihood functionの形態をみる

logL <- function(m) sum(dpois(data, m, log = TRUE))
lambda <- seq(2, 5, 0.1)
plot(lambda, sapply(lambda, logL), type = "b")

L <- function(m) prod(dpois(data, m, log = FALSE))
plot(lambda, sapply(lambda, L), type = "b")

最尤法を用いるポアソン回帰を行って切片のexp()を求めると同じ値

fit.poisson <- glm(plant.seed ~ 1, family = poisson)
exp(coef(fit.poisson))

(Intercept) 
       3.56