Rで「フィールドデータによる統計モデリングとAIC」のP.95のシュミレーションをしてみた

鳴らない携帯に月7000円を払うなど.

あっ, どーも僕です.

現実が地獄でとにかく逃避したくてタイトルにあるシュミレーションをRでやってみました.

フィールドデータによる統計モデリングとAIC (ISMシリーズ:進化する統計数理)

シュミレーションの目的は「最大対数尤度と平均対数尤度の差を計算し, その平均値がパラメータ数で近似できるということを確認する」です.

このことが確認できれば, AICの理解に役立つということですね.

このシュミレーションは真のモデルがわかっているという前提で, 多項式近似の正規分布モデルを当てはめていこうというのが流れです.

実際には真のモデルはわからないのであくまでシュミレーションですね.

で, 具体的になにをするのかというと.

まず真のモデルは正規分布モデルであり, 説明変数xに対して目的変数yが次の式を平均値として標準偏差3でばらつくと考えます.

真のモデルのイメージはこちら. ピンクが平均値で破線が1.96σです.

でもって, 観測したデータから次の式から得る最大対数尤度

と次の式から得る平均対数尤度の差を計算します.

このとき, 観測したデータに当てはめるモデルは次に示すような多項式正規分布モデルで1次から6次でやります.

まっ, 詳しく説明するのはしんどいので本みてください・・・・すんません.

実際にRでシュミレーションをしたスクリプトがこちら. ほんと同様にサンプル数200でやっています. このシュミレーションはループを回さなくても行列演算で一度に計算できます.シュミレーション回数は1000回です.

# ***************************************************
# 最大対数尤度と平均対数尤度の差のシュミレーション
# **************************************************
# 真のモデル4次関数に正規分布の誤差をもたせる
# 真のモデルはわかっている状態で正規分布モデルをつくる
# 多項式近似で1次から6次まで最尤推定量と平均対数尤度をの差を計算



#true.model
true.model <- function(x) x^4 - 2 * x^3 - 3 * x^2

#データ： サンプル数200 シュミレーション回数1000
#このスクリプトだとN=1000くらいが限界
N <- 200
N.sim <- 1000
sigma_0 <- 3
x <- seq(0.8, 2, length = N)
Y <- sapply(1:N.sim, function(i, x) true.model(x) + rnorm(n = N, m = 0, sd = sigma_0 ), x)

#デザインマトリックスだっけ？
X <- matrix(c(rep(1, N), x, x^2 , x^3, x^4, x^5, x^6), nrow = N)

#係数の最尤推定量
Estimate_a <- 
  lapply(1:6, function(i) solve(crossprod(X[,1:(i+1)]), crossprod(X[,1:(i+1)], Y)))

#分散の最尤推定量
Estimate_sigma <- 
  sapply(1:6, function(i) colMeans((Y - X[,1:(i+1)] %*% Estimate_a[[i]])^2))

#最大対数尤度
MaximumLogLikelihood <- - N * ( 1 + log(2 * pi * Estimate_sigma)) / 2

#平均対数尤度
error <- 
  sapply(1:6, function(i) colSums((true.model(x)- X[,1:(i+1)] %*% Estimate_a[[i]])^2))
MeanLogLikelihood <-
  -(N * log(2 * pi * Estimate_sigma) + N * sigma_0^2 / Estimate_sigma + error / Estimate_sigma) * 0.5

#最大対数尤度と平均対数尤度の差を求める
result <- MaximumLogLikelihood - MeanLogLikelihood

結果をグラフ化したものが次です. パラメータ数は多項式の次数+定数項+分散, つまり次数＋２です. 確かに最大対数尤度と平均対数尤度の差の平均値が, パラメータ数で近似できそうなことがわかります.

最大対数尤度と平均対数尤度の差の分布はこちら. 平均対数尤度が非対称な分布なのでこの分布も非対称ですね.

最後に, グラフづくりのスクリプトを.

#作図
#パラメータ数と平均差
def.par <- par(no.readonly = TRUE) # デフォルトの環境を保存
par(mai=c(1.5,1.5,1,1))
plot(0, xlim = c(2,9), ylim = c(2,9), axes = F, type = "n", yaxs ="i",xaxs="i",
     xlab = "Number of parameters", ylab = "average difference", 
     cex.lab = 1.2, font.lab = 4)
usr <- par("usr")
rect(usr[1], usr[3], usr[2], usr[4], col = "#F8F8F8")
abline(h=0:9, v=0:9, col = "lightgrey")
axis(side = 1, tcl = 0.4, at = 2:9, cex.axis = 1.2)
axis(side = 2, tcl = 0.4, at = 2:9, cex.axis = 1.2)
abline(a = 0, b = 1, lty = 5)
box()
points(3:8, colMeans(result), col = "#0068b7", type ="b", cex = 2, pch = 19)
par(def.par)　　　　　　　　　　　 # デフォルトの環境を復帰



#差の分布
title=paste("Following regression N = ", 1:6, sep = "")
layout(matrix(1:6,ncol=3))
myplot <- function(i, X) plot(density(X[,i]), lwd = 3, type ="l", tcl = 0.4, 
                              cex.axis = 1.2, main = title[i], col = "#0068b7", 
                              xlab = "", ylab = "")
sapply(1:6, myplot, result)