モデルのパラメータの分布を眺めてみる

今日の時系列の授業は時系列というより統計っぽい話。最小二乗法とかからパラメータの分布は自分で一回見てみないとあれだよ、って言われたので調べといた。

説明変数2つくらいのモデルを考える。10000個のデータセットを100個ずつの100個のデータに分ける。それぞれのデータから一つモデルを作り、モデルを計100個作る。それぞれのモデルのパラメータを引っ張ってきて密度トレイスを描いてみた。正規…分布?


適当すぎる実験コードは下。方法がそもそもこれであっているのがが危ない。てか、for文が遅いのでapplyファミリー走らせないと。。。

#計量時系列分析のやつで遊ぶ

#二変数くらいでモデルを考えて、適当に誤差項でばらまく

y <- c()

n <- 10000

x1 <- c()
x2 <- c()
rnorm(1)

for(i in 1:n){
  x1[i] <- rnorm(1)*10
  x2[i] <- rnorm(1)*3
  y[i] <- 3*x1[i]+10*x2[i]+rnorm(1)*5
}
summary(y)

x <- cbind(x1,x2)
summary(x)

beta <- matrix(0,2,100)

beta

#とりあえずこれで出る
solve(t(x)%*%x)%*%t(x)%*%y

seq(1,10000,by=100)

for(i in seq(1,10000,by=100)){
  print(i+1)
}
9901+99
y[10000+1]

warnings()
count <- 1
for(i in seq(1,10000,by=100)){
  beta[,count] <- solve(t(x[i:(i+99),])%*%x[i:(i+99),])%*%t(x[i:(i+99),])%*%(y[i:(i+99)])
  count <- count+1
}
beta

hist(beta[1,],nclass=20)
hist(beta[2,],nclass=20)
png("den1.jpg")
plot(density(beta[1,]))
dev.off()
png("den2.jpg")
plot(density(beta[2,]))
dev.off()
mean(beta[1,])
mean(beta[2,])

var(beta[1,])
var(beta[2,])

25*solve(t(x)%*%x)