emp_alpha<-function(nsims){
n<-20 # sample size of 20
tc<-qt(1-0.05,n-1) #this is 1.729, i.e the reject region,here we suppose u=0
ic<-0
for(i in 1:nsims){
samp<-rlogis(n,0,1) # generate 20 observations from logistic distribution with u=0
ttest<-(sqrt(n)*mean(samp))/var(samp)^.5 # calcuate ts from sample.
if(ttest > tc){ic<-ic+1}
}
alpha<-ic/nsims
err<-1.96*sqrt((alpha*(1-alpha))/nsims)
list(empiricalalpha=alpha,error=err)
}
emp_alpha(100000)
#######simulate a power function
emp_power<-function(nsims,u){
n<-20 # sample size of 20
tc<-qt(1-0.05,n-1) #this is 1.729, i.e the reject region,here we suppose u=0
ic<-0
for(i in 1:nsims){
samp<-rlogis(n,u,1) # generate 20 samples from the logistic distribuiton with u
ttest<-(sqrt(n)*mean(samp))/var(samp)^.5 # calcuate ts from sample.
if(ttest > tc){ic<-ic+1}
}
power<-ic/nsims
return(power)
#err<-1.96*sqrt((power*(1-power))/nsims)
#list(emppower=power,error=err)
}
emp_power(1000,0.5) #power for u=0.5
#We will draw a power function for this simulation
powers<-NULL
u<-seq(0,1,0.001)
for(i in 1:1001){
powers[i]<-emp_power(1000,u[i])
}
plot(u,powers,type="l") # this is power function plot
