**MuMInパッケージを利用するので、事前にインストールが必要
真のモデル:
yは平均μ、σ=sdの正規分布に従い
μ=1.5X1-1.1X2+20
*簡単のため、X1とX2は直交しているものとする。
データ：
真のモデルからn個のサンプルをランダムに抽出
推定：
線形回帰モデルを作り、n個のサンプルを使ってパラメータを最尤推定する
変数選択：
y~X1+X2+Intercept
をフルモデルと見なし、情報量基準が最小のモデルを”ベストモデル”として採用する
結果の格納：
ベストモデルにおける各パラメータの最尤推定値を保存する
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
以上のプロセスを1万回反復することで、post-model selectionにおける最尤推定値の挙動を確認できる。
通常、最尤推定量が”良い”推定量とみなされるのは
1)漸近正規性
2)普遍性
3)一致性
という三つの性質を備えていることが根拠であるが、モデル選択後の推定量では1や2の性質が容易に失われることが確認出来る。
##関数の作成
BIC_sim<-function(n,sd){
library(MuMIn)
options(na.action = "na.fail")
x1<-matrix(0,nrow=10000,ncol=n)
x2<-matrix(0,nrow=10000,ncol=n)
y<-matrix(0,nrow=10000,ncol=n)
for (i in 1:10000){
x1[i,]<-runif(n,0,100)
x2[i,]<-rnorm(n,0,10)
y[i,]<-1.5*x1[i,]-1.1*x2[i,]+20+rnorm(n,mean=0,sd=sd)
}
temp_x<-matrix(0,nrow=10000,ncol=3)
for (i in 1:10000){
temp1<-lm(y[i,]~x1[i,]+x2[i,])
temp2<-dredge(temp1,rank="BIC",trace=FALSE) ###引数rankを変更することで、他の情報量基準における挙動も確認できる
temp_x[i,1]<-temp2[1,1]
temp_x[i,2]<-temp2[1,2]
temp_x[i,3]<-temp2[1,3]
}
return (temp_x)
}
##シミュレーションの実行
a<-(BIC_sim(100,300))
##最尤推定値の確認
summary(a)
hist(a[,2])