回帰解析ジョブ1

1620 ワード

作業内容:
1、『線形統計モデル』(王松桂等)の再現例3.1.3
2、シミュレーションの最小二乗法、そしてシミュレーションの結果に基づいて推定係数と推定分散のいくつかの結論を出して、例えば偏りがないなど
ジョブインプリメンテーション(リファレンスのみ)
#1
y<-c(10.98,11.13,12.51,8.40,9.27,8.73,6.63,8.50,7.82,9.14,8.24,12.19,11.8,
9.57,10.94,9.58,10.09,8.11,6.83,8.88,7.68,8.47,8.86,10.36,11.08)
x<-c(35.30,29.70,30.80,58.8,61.4,71.3,74.4,76.7,70.7,57.5,46.4,28.90,28.1,
39.1,46.80,48.5,59.30,70.0,70.0,74.5,72.1,58.1,44.6,33.40,28.60)
lm(y~x)
summary(lm(y~x))
Rss<-deviance(lm(y~x))
Rss/23
#         :
#beta0=13.57290,beta1=-0.07873,sigma^2=0.760774
#     
x1<-x-mean(x)
lm(y~x1)
summary(lm(y~x1))
Rss<-deviance(lm(y~x1))
Rss/23
#         :
#beta0= 9.43160,beta1=--0.07873,sigma^2=0.760774
plot(y~x)
abline(lm(y~x))#see picture1
plot(y~x1)
abline(lm(y~x1))#see picture2


#2    :y=3.2*x+0.5+e
x<-rexp(100,3)
simulation<-function(x){
beta0<-rep(0,1000)
beta1<-rep(0,1000)
sigma<-rep(0,1000)
Rss<-rep(0,1000)
for(i in 1:1000){
set.seed(7*i)
e<-rnorm(100)
y<-3.2*x+0.5+e
beta0[i]<-lm(y~x)$coef[1]
beta1[i]<-lm(y~x)$coef[2]
Rss[i]<-deviance(lm(y~x))
sigma[i]<-deviance(lm(y~x))/98
}
print("beta0=")
print(mean(beta0))
print("beta1=")
print(mean(beta1))
print("sigma^2=")
print(mean(sigma))
ks.test(Rss,"pchisq",98)
}
simulation(x)
#    :
#[1] "beta0="
#[1] 0.4965268
#[1] "beta1="
#[1] 3.199646
#[1] "sigma^2="
#[1] 1.003083
#
#        One-sample Kolmogorov-Smirnov test
#data:  Rss 
#D = 0.0231, p-value = 0.6619
#alternative hypothesis: two-sided 
#    :     3.2.1,3.2.3