gofsim=function(y,X,intc,msim){
if(intc==1)lm1=lm(y~X)
if(intc==0)lm1=lm(y~X-1)
#r1=rstandard(lm1)
r1=lm1$resid
n=length(r1)
q1<-qqnorm(r1,plot=F)
c1<-cor(q1$x,q1$y)
u1<-pnorm(sort(r1),mean=0,sd=summary(lm1)$sigma)
d1<-max(c((1:n)/n-u1,u1-(0:(n-1))/n))
w1<-sum((u1-((1:n)-0.5)/n)^2)+1/(12*n)
a1<--sum((2*(1:n)-1)*log(u1)+(2*n+1-2*(1:n))*log(1-u1))/n-n
count<-rep(0,4)
for(j in 1:msim){
ysim<-rnorm(n)
if(intc==1)lm2=lm(ysim~X)
if(intc==0)lm2=lm(ysim~X-1)
#r2=rstandard(lm2)
r2=lm2$resid
q2<-qqnorm(r2,plot=F)
c2<-cor(q2$x,q2$y)
u2<-pnorm(sort(r2),mean=0,sd=summary(lm2)$sigma)
d2<-max(c((1:n)/n-u2,u2-(0:(n-1))/n))
w2<-sum((u2-((1:n)-0.5)/n)^2)+1/(12*n)
a2<--sum((2*(1:n)-1)*log(u2)+(2*n+1-2*(1:n))*log(1-u2))/n-n
if(c2<c1)count[1]<-count[1]+1
if(d2>d1)count[2]<-count[2]+1
if(w2>w1)count[3]<-count[3]+1
if(a2>a1)count[4]<-count[4]+1
}
ddf=matrix(c(c1,d1,w1,a1,count/msim),byrow=T,ncol=4)
ddf=round(ddf,5)
ddf=cbind(c('Test','p-value'),ddf)
ddf=rbind(c('','LG','KS','CVM','AD'),ddf)
ddf
}