Friday, September 23, 2005

using R : regression test - 회귀분석에 대한 모의실험

y=a+bx+e
x=1:n
y=alpha+beta*x+rnorm(n)
>lm.out=lm(y~x)
#lm.out 은 lm()에 의해 생성된 객체, 적합된 결과는 객체인 lm.out을 summary()함수의 매개변수로 지정함으로써 확인가능
>summary(lm.out)
#회귀계수들의 추정값은 lm.out$coefficients
>lm.out$coefficients
#추정된 회귀계수들의 분산값은 vcov()함수를 이용
>var=vcov(lm.out)
>var
>var[1,1]

#Reg.confidence=function(nt,n,beta)
#beta = vector of true regression coeff
#nt=No.of replication
#n=sample size
{win.graph()
par(mfrow=c(2,1))
mtitle=c("Intercept","Slope")
x=1:n
Exb=beta[1]+beta[2]*x
t.val=qt(0.975,n-2)
trial.a=trial.b=rep(0,nt)

for(i in 1:nt)
{ y=Exb+rnorm(n)
lm.out=lm(y~x)
std=sqrt(diag(vcov(lm.out))
if(abs(lm.out$coefficients[1]-beta[1])if(abs(lm.out$coefficients[2]-beta[2])}
trial.a=cumsum(trial.a)/1:nt
plot(trial.a,type="l",main=mtitle[1],ylim=c(0.9,1))
abline(0.95,0)
trial.b=cumsum(trial.b)/1:nt
plot(trial.b,type="l",main=mtitle[2],ylim=c(0.9,1))
abline(0.95,0)
}

No comments: