SSブログ

スタイン推定量を計算してみた [統計]

『ベイズモデリングの世界』講義1にでてくるスタイン推定量を計算してみて、誤差平方和の期待値がちいさくなることを確認してみたメモです。

Rコードです。まず、残差平方和とスタイン推定量の関数を定義しています。

## Sum of squared error
sse <- function(theta.star, theta) {
  sum((theta.star - theta)^2)
}

## Stein estimator
stein <- function(y, sigma) {
  n <- length(y)
  if (n <= 3)
    stop("n must be larger than 3.")
  s2 <- sum((y - mean(y))^2) / (n - 3)
  a <- sigma^2 / s2
  theta.stein <- (1 - a) * y + a * mean(y)
}

まず、θを[0, 128]の一様乱数で生成させてみました。σは2としました。

set.seed(123)
n <- 64
theta <- runif(n, 0, 128)
sigma <- 2
sse1 <- sapply(1:2000, function(i) {
  y <- rnorm(n, theta, sigma)
  theta.stein <- stein(y, sigma)
  c(sse(y, theta), sse(theta.stein, theta))
})
apply(sse1, 1, mean)

結果です。たしかにスタイン推定量の方が誤差平方和がちいさくなりました。

[1] 255.7266 254.9792

つぎに、θをNormal(0, 1)の正規乱数で生成させてみました。σはやはり2です。

set.seed(123)
n <- 64
theta <- rnorm(n, 0, 1)
sigma <- 2
sse2 <- sapply(1:2000, function(i) {
  y <- rnorm(n, theta, sigma)
  theta.stein <- stein(y, sigma)
  c(sse(y, theta), sse(theta.stein, theta))
})
apply(sse2, 1, mean)

結果です。スタイン推定量を使った方がかなりちいさくなりました。「関係ある」データほど、平均の方に縮小するということです。

[1] 255.74957  52.75787

最後に、ほぼ「関係のない」データとなるようにしてみました。

set.seed(123)
n <- 64
theta <- runif(n, -1e+5, 1e+5)
sigma <- 1
sse3 <- sapply(1:2000, function(i) {
  y <- rnorm(n, theta, sigma)
  theta.stein <- stein(y, sigma)
  c(sse(y, theta), sse(theta.stein, theta))
})
apply(sse3, 1, mean)

この場合、誤差平方和の平均に差はほぼありませんでした。

[1] 63.93165 63.93165

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

nice! 3

コメント 0

コメントを書く

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

Facebook コメント