So-net無料ブログ作成

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

分散の不変推定量とか、最尤推定値の漸近正規性とか。

分散の不変推定量はN-1でわることを、シミュレーションで確認する。

### 2.8
## 分散の不変推定量
set.seed(1)
N <- 20
lambda <- 4
R <- 1000
x <- replicate(R, {
  x <- rpois(N, lambda)
  v1 <- sum((x - mean(x))^2) / N
  v2 <- sum((x - mean(x))^2) / (N - 1)
  c(v1, v2)
})
apply(x, 1, mean)
[1] 3.801732 4.001824

最尤推定値は、サンプルサイズがおおきくなるにつれ、真値に集中してくる。

### 2.10
sim2 <- function(R = 100, N = c(4, 16), lambda = 4, binwidth = NULL) {
  labels <- sapply(N, function(n) paste("n =", n))
  x <- matrix(0, length(N), R)
  for (i in seq_along(N)) {
    x[i, ] <- replicate(R, {
      ## Maximum Likelihood Estimate
      d <- rpois(N[i], lambda)
      optim(1,
            fn = function(lambda) {
              lp <- sum(log(dpois(d, lambda)));
              return(-lp)
            },
            method = "Brent",
            lower = 1.0e-4, upper = 1.0e+3)$par
    })
  }
  data.frame(label = factor(rep(labels, R),
                            levels = labels),
             x = c(x)) %>%
    ggplot() +
    geom_freqpoly(aes(x = x, color = label), binwidth = binwidth) +
    labs(x = "最尤推定値", y = "回数") +
    scale_y_continuous(expand = c(0, 0)) +
    scale_color_discrete(name = "") +
    theme_classic(base_family = "IPAexGothic") +
    theme(legend.justification = c(1, 1), legend.position = c(1, 1))
}

sim2(R = 1000, N = c(4, 16, 64, 256), binwidth = 0.4)
sim2(R = 1000, N = c(256, 1024, 4096, 16384), binwidth = 0.05)

2.5a.png
2.5b.png

最尤推定値は、サンプルサイズがおおきいと、正規分布にちかづく。

### 2.11
sim3 <- function(R = 100, N = 4, lambda = 4, binwidth = 0.4) {
  label <- paste("n =", N)
  x <- replicate(R, {
    ## Maximum Likelihood Estimate
    d <- rpois(N, lambda)
    optim(1,
          fn = function(lambda) {
            lp <- sum(log(dpois(d, lambda)));
            return(-lp)
          },
          method = "Brent",
          lower = 1.0e-4, upper = 1.0e+3)$par
  })
  ## 漸近正規推定量
  mu <- lambda
  sigma <- sqrt(lambda / N)
  fun <- function(x, R, mu, sigma, binwidth) {
    R * binwidth * dnorm(x, mu, sd = sigma)
  }
  data.frame(x = x) %>%
    ggplot() +
    geom_histogram(aes(x = x), alpha = 0.5, binwidth = binwidth) +
    stat_function(fun = fun,
                  args = list(R = R, mu = mu, sigma = sigma, binwidth = binwidth)) +
    labs(title = label, x = "最尤推定値", y = "回数") +
    scale_y_continuous(expand = c(0, 0)) +
    theme_classic(base_family = "IPAexGothic")
}

set.seed(1)
sim3(R = 10000, N = 4, binwidth = 0.4)
sim3(R = 10000, N = 64, binwidth = 0.2)
sim3(R = 10000, N = 1024, binwidth = 0.05)
sim3(R = 10000, N = 16384, binwidth = 0.01)

2.6a.png
2.6b.png
2.6c.png
2.6d.png


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

nice! 1

コメント 0

コメントを書く

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

Facebook コメント