#determine the threshold for the test H0: p=2/3 versus H1: p=1/3 plot(0:100,pbinom(0:100,size=100,prob=2/3),type="p") abline(h=0.05) t<-qbinom(p=0.05,size=100,prob=2/3)-1 pbinom(t,size = 100,prob=2/3) #calculate beta 1-pbinom(t,size=100,prob=1/3) #evaluate performance of the developed test # one run under H0: n<-100 p<-2/3 x<-sample(c(1,0),size=n,replace = TRUE,prob=c(2/3,1/3)) if(length(x[x==1])<=58){print("reject H0")} # R=10000 runs under H0 R<-10000 reject<-rep(0,R) for(i in 1:R){ x<-sample(c(1,0),size=n,replace = TRUE,prob=c(2/3,1/3)) if(length(x[x==1])<=58){reject[i]<-1} } mean(reject) barplot(table(reject)) # R=10000 runs under H1 #-------------------------------------------------------------------- #one sample t-test manually; test mu=mu0 versus two-sided alternative n <- 20 mu0 <- 3 sigma <- 1 x <- rnorm(20,mean=mu0,sd=sigma) alpha <- 0.05 Tn <- sqrt(n)*(mean(x)-mu0)/sd(x) if(abs(Tn)> qt(1-alpha/2,df=n-1)){print("reject H0")} #systematic - check error of first type n <- 20 mu0 <- 3 sigma <- 1 alpha <- 0.05 R <- 1000 errorI <- rep(0,R) for(i in 1:R){ x <- rnorm(n,mean=mu0,sd=sigma) Tn <- sqrt(n)*(mean(x)-mu0)/sd(x) if(abs(Tn)> qt(1-alpha/2,df=n-1)){errorI[i]<-1} } mean(errorI) #check power of the test for the case mu=2 n <- 20000 mu0 <- 3 mu <- 2.9 sigma <- 1 alpha <- 0.05 R <- 1000 errorII <- rep(0,R) for(i in 1:R){ x <- rnorm(n,mean=mu,sd=sigma) Tn <- sqrt(n)*(mean(x)-mu0)/sd(x) if(abs(Tn)<= qt(1-alpha/2,df=n-1)){errorII[i]<-1} } mean(errorII) power <- 1- mean(errorII) power