So-net無料ブログ作成

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

カルバック・ライブラー情報量とか、平均対数尤度と最大対数尤度との差とか。

強度の異なるポアソン分布を比較する。

library(ggplot2)
library(dplyr)

### 4.2
k <- 0:20
p <- dpois(k, 4)
q1 <- dpois(k, 3)
q2 <- dpois(k, 5)
q3 <- dpois(k, 3.3)
q4 <- dpois(k, 4.8)

## p, q1, q2
data.frame(k = rep(k, 3),
           x = c(p, q1, q2),
           l = rep(c("強度4", "強度3", "強度5"),
                   each = length(k))) %>%
  ggplot() +
  geom_line(mapping = aes(x = k, y = x, colour = l)) +
  scale_color_discrete(name = "") +
  theme_classic(base_family = "IPAexGothic") +
  theme(legend.justification = c(1, 1), legend.position = c(1, 1))

4.1a.png

data.frame(k = rep(k, 3),
           x = c(p, q3, q4),
           l = rep(c("強度4", "強度3.3", "強度4.8"),
                   each = length(k))) %>%
  ggplot() +
  geom_line(mapping = aes(x = k, y = x, colour = l)) +
  scale_color_discrete(name = "") +
  theme_classic(base_family = "IPAexGothic") +
  theme(legend.justification = c(1, 1), legend.position = c(1, 1))

4.1b.png

カルバック・ライブラー情報量を定義

## Kullback-Leibler information
KL <- function(p, q) {
  sum(p * (log(p) - log(q)))
}

カルバック・ライブラー情報量、差の平方和、差の絶対値の和

> KL(p, q1)
[1] 0.1507283
> KL(p, q2)
[1] 0.1074258
> KL(p, q3)
[1] 0.06948756
> KL(p, q4)
[1] 0.07071378
> 
> sum((p - q1)^2)
[1] 0.02237314
> sum((p - q2)^2)
[1] 0.01515742
> sum((p - q3)^2)
[1] 0.01038763
> sum((p - q4)^2)
[1] 0.01011685
> 
> sum(abs(p - q1))
[1] 0.4275235
> sum(abs(p - q2))
[1] 0.3766872
> sum(abs(p - q3))
[1] 0.2937362
> sum(abs(p - q4))
[1] 0.3051563

差の平方和と差の絶対値の和の値が本とは違っていたりするが。

平均対数尤度と最大対数尤度との差が、パラメーター数に一致することをシミュレーションで確認する。

### 4.3
set.seed(1)
n <- 100
lambda <- 4
k <- 0:50

## mean log-likelihood
mean.loglik <- function(n, p, q) {
  n * sum(p * log(q))
}

sim1 <- function(R = 100, n = 100, lambda = 4,
                 k = 0:50) {
  p <- dpois(k, lambda)
  z <- replicate(R, {
    x <- rpois(n, lambda)
    ml <- optim(1,
                fn = function(lambda) {
                  -sum(log(dpois(x, lambda)))
                },
                method = "Brent",
                lower = 0.0001, upper = 1000)
    l <- -ml$value
    r <- ml$par
    q <- dpois(k, r)
    m <- mean.loglik(n, p, q)

    l - m
  })
  return(z)
}

mean(sim1(R = 100000))
[1] 0.9965566

だいたい1になった。

パラメーター数が2のとき。

### 4.5

sim2 <- function(R = 100, n1 = 50, n2 = 50,
                 lambda1 = 4, lambda2 = 6,
                 k = 0:50) {
  p1 <- dpois(k, lambda1)
  p2 <- dpois(k, lambda2)
  z <- replicate(R, {
    x <- c(rpois(n1, lambda1), rpois(n2, lambda2))
    ml <- optim(c(2, 2),
                fn = function(lambda) {
                  -sum(log(dpois(x[1:n1], lambda[1]))) -
                    sum(log(dpois(x[(n1 + 1):(n1 + n2)], lambda[2])))
                })
    l <- -ml$value
    r <- ml$par
    q1 <- dpois(k, r[1])
    q2 <- dpois(k, r[2])
    m <- mean.loglik(n1, p1, q1) + mean.loglik(n2, p2, q2)

    l - m
  })
  return(z)
}

mean(sim2(R = 10000))
[1] 1.958838

だいたい2になった。


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

nice! 2

コメント 0

コメントを書く

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

Facebook コメント