So-net無料ブログ作成

『ポアソン分布・ポアソン回帰・ポアソン過程』をRで再現してみる(3) [統計]

グラフにかぎらず、再現してみようかと。

尤度を計算する。

library(ggplot2)
library(dplyr)

### 2.3
## 本数
N <- 15
## 面積
A <- 143 * 139 / 10000

## 尤度関数
lik.fun <- function(N, r, A) {
  dpois(N, r * A)
}

## 密度と尤度のグラフ
ggplot(data = data.frame(r = 0), mapping = aes(x = r)) +
  stat_function(fun = lik.fun, args = list(N = N, A = A)) +
  xlim(5, 10) +
  ylim(0, 0.12) +
  labs(x = "密度r",
       y = "尤度L(r)") +
  theme_classic(base_family = "IPAexGothic")

2.2.png

最尤推定する。

## 最尤推定
ml <- optim(1, fn = function(r, N, A) -dpois(N, r * A),
            lower = 0, upper = 100,
            method = "Brent",
            N = N, A = A)
r <- ml$par
print(r)
[1] 7.54641

尤度と対数尤度

### 2.4
rm(list=ls())
library(ggplot2)
library(dplyr)

## ホオノキの本数
Mo <- c(1, 2, 3, 3, 1, 0,
        0, 2, 0, 0, 0, 2,
        0, 0, 3, 2, 0, 1)

## 尤度関数
lik.fun <- function(N, r, A) {
  l <- 1
  for (i in 1:length(N))
    l <- l * dpois(N[i], r * A)
  return(l)
}

ggplot(data = data.frame(r = 0), mapping = aes(x = r)) +
  stat_function(fun = lik.fun, args = list(N = Mo, A = 0.5)) +
  xlim(0, 6) +
  labs(x = "密度r",
       y = "尤度L(r)") +
  theme_classic(base_family = "IPAexGothic")

2.3a.png

## 対数尤度関数
loglik.fun <- function(N, r, A) {
  loglik <- 0
  for (i in seq_along(N)) {
    loglik <- loglik + log(dpois(N[i], r * A))
  }
  return(loglik)
}

ggplot(data = data.frame(r = 0), mapping = aes(x = r)) +
  stat_function(fun = loglik.fun, args = list(N = Mo, A = 0.5)) +
  xlim(0.01, 6) +
  labs(x = "密度r",
       y = "尤度L(r)") +
  theme_classic(base_family = "IPAexGothic")

2.3b.png

最尤推定

## 最尤推定
ml2 <- optim(1, fn = function(N, r, A) -loglik.fun(N, r, A),
             method = "Brent",
             lower = 0.001, upper = 100,
             N = Mo, A = 0.5)
ml2$par
[1] 2.222222

95%信頼区間にデータがおさまっているか。

### 2.5
set.seed(1)

sim1 <- function(data, R, n.max = max(data)) {
  N <- length(data)

  ## Maximum Likelihood Estimate
  r <- optim(1,
             fn = function(lambda) {
               lp <- sum(log(dpois(data, lambda)));
               return(-lp)
             },
             method = "Brent",
             lower = 0.001, upper = 100)$par
  x <- sapply(1:R, function(i) {
    n <- rpois(N, r)
    table(factor(n, levels = 0:n.max))
  })
  c <- as.numeric(table(factor(data, levels = 0:n.max)))
  m <- apply(x, 1, mean)
  u <- apply(x, 1, quantile, probs = 0.975)
  l <- apply(x, 1, quantile, probs = 0.025)
  data.frame(n = 0:n.max,
             count = c, mean = m, upper = u, lower = l) %>%
    ggplot() +
    geom_ribbon(aes(x = n, ymin = lower, ymax = upper), alpha = 0.33) +
    geom_line(aes(x = n, y = mean), linetype = 2) +
    geom_line(aes(x = n, y = count)) +
    geom_point(aes(x = n, y = count)) +
    ylim(range(pretty(0:max(u)))) +
    labs(x = "本数k", y = "調査区の数") +
    theme_classic(base_family = "IPAexGothic")
}
sim1(Mo, 1000, 8)

2.4a.png

## ナナカマドの本数
Sc <- c(3, 3, 2, 3, 7, 10,
        0, 0, 0, 1, 1, 2,
        0, 1, 1, 1, 4, 0)
sim1(Sc, 1000, 11)

2.4b.png


タグ:R
nice!(1)  コメント(0) 
共通テーマ:日記・雑感

nice! 1

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

Facebook コメント