スタイン推定量を計算してみた [統計]
『ベイズモデリングの世界』講義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
コメント 0