rm(list=ls())
N=25
b0=1
b1=1.5
trials=5000
bst=vector(,trials)
se=vector(,trials)
tstore=vector(,trials)
reject=vector(,trials)

for (i in 1:trials){
	x1=rnorm(N)
	y=b0+b1*x1+rnorm(N)
	output=lm(y~x1)
	m=summary(output)
	b2=m$coef
      yhat=m$fitted
      uhats=m$residuals
	bst[i]<-b2[2,1]
	se[i]<-b2[2,2]
	tstore[i]=(bst[i]-b1)/se[i]
      

	#crit=qt(.025,N-2)
      crit=1.96
	reject[i]=(abs(tstore[i])>=abs(crit))
      
}
hist(bst)
print("Rejection Probability Under H0")
sum(reject)/trials