As we know, the power of a test is the probability of given H1 being true, and rejecting H0.
Usually, we do not know exactly the power of a test, because we are not sure where the true value of the estimated parameter is. But for some specific situations, we can tell how much the power is, just like this example:
Assume we have n i.i.d. samples of X ~ N(mu,sigma^2), and alpha = 0.05.
H0: mu = 0
H1: mu = 0.01, 0.02, ...... 0.98, 0.99, 1.00
We have 3 kinds of n, n = 15, 30, 50
In the Monte Carlo simulation, we generate X with normal distribution, which has mu corresponding to H1, and sigma = 1.
Then we construct Z statistics which follows standard normal distribution with given sigma :
Z = n^(1/2)*(X_bar - 0)/sigma
Here we reject H0 if Z >z0.95, thus, we are doing a one-sided test, only the values which are abnormally large will be rejected.
With such setting, if we let D = 1 when rejecting H0, and D = 0 if not reject, for each mu and n, the power function Beta = mean of D could represent the power of this test. (with 10000 times simulations)
Note that, because here we use z0.95 as critical value, we are trying to find the tail probability of a series of normal distributions given critical value z0.95 ( this z0.95 is for standard normal distribution).
In addition, if we assume we do not know the sigma, and we have to estimate it with sample standard deviation estimator, then the test statistic Z follows t distribution, but here we still compare it with normal distribution results.
In the figure above, the dashed lines are theoretical power functions, the lines are the simulated power of the test, and the crosses are the simulated power of the test using estimated sigma.
It is not hard to find that the performances of Z given a known sigma are good, and Z with estimated sigma do not have such a good fit.
rm(list=ls())
R <- 10000
mu.X.true.seq <- seq(0,1,len=100)
beta.n <- list(matrix(NA,100,R),matrix(NA,100,R),matrix(NA,100,R))
z_crit <- qnorm(0.95,mean=0,sd=1)
N <- c(15,30,50)
for (m in 1:3){
for (i in 1:100) {
for (j in 1:R) {
X <- rnorm(N[m],mu.X.true.seq[i],1)
Judge <- if(sqrt(N[m])*mean(X)>z_crit){
1}else{0}
beta.n[[m]][i,j]<- Judge
}
}}
result <- matrix(NA,100,3)
for (i in 1:3) {
for (j in 1:100) {
result[j,i] <- mean(beta.n[[i]][j,])
}
}
Gauss.beta<- function(n, ###sample size
alpha=0.05, ####significance level
mu.X.true=0, ###true mean of the X_i
mu.X.null=0, ###mean under the null hypothesis
var.X =1 ###assumed (known) variance of X_i
){
###Critical value
z_crit <- qnorm(1-alpha,mean=0,sd=1)
#power
Phi<-pnorm(q =z_crit, mean =sqrt(n)*(mu.X.true-mu.X.null)/sqrt(var.X), sd=1)
power <-1-Phi
#Return result
return(power)
}
##Vectorisation
mu.X.true.seq <- seq(0,1,len=100)
## Trajectories of the power function for different sample sizes:
## n=15
beta.n.15<- Gauss.beta(n=15, mu.X.true=mu.X.true.seq)
## n=30
beta.n.30<- Gauss.beta(n=30, mu.X.true=mu.X.true.seq)
## n=50
beta.n.50<- Gauss.beta(n=50, mu.X.true=mu.X.true.seq)
###plot
par(mar=c(5.1,4.1+1,4.1,2.1))
plot(y=0, x=0, type= "n",
ylim = c(0,1),
xlim=range(mu.X.true.seq),
xlab =expression(paste(mu," (True Mean of ",X[i],")")),
ylab = expression(paste(" Power function ",beta[n]^Z,(mu))),
main = "Power function of the (one-sided) Gauss test")
##Null hypothesis mean
mtext(text = expression(mu[0]==0), side =1, line =2, at = 0)
##Trajectories
lines(y=beta.n.15, x=mu.X.true.seq, lty =2,lwd=2)
lines(y=beta.n.30, x=mu.X.true.seq, lty =3,lwd=2)
lines(y=beta.n.50, x=mu.X.true.seq, lty =4,lwd=2)
lines(mu.X.true.seq,result[,1],lwd=1.5)
lines(mu.X.true.seq,result[,2],col="red",lwd=1.5)
lines(mu.X.true.seq,result[,3],col="blue",lwd=1.5)
###significance level
axis(4, at =0.05, labels= expression(alpha == 0.05))
abline (h=0.05, lty =1)
###Legend
legend("topleft", title = "Sample Sizes", legend = c("n=50", "n=30", "n=15"), lty = c(4:2)
,lwd=rep(2,3))
############################################
mu.X.true.seq <- seq(0,1,len=100)
beta.n2 <- list(matrix(NA,100,R),matrix(NA,100,R),matrix(NA,100,R))
z_crit <- qnorm(0.95,mean=0,sd=1)
N <- c(15,30,50)
for (m in 1:3){
for (i in 1:100) {
for (j in 1:R) {
X <- rnorm(N[m],mu.X.true.seq[i],1)
Judge <- if(sqrt(N[m])*mean(X)/(1/(N[m]-1)*sum((X-mean(X))^2))^0.5>z_crit){
1}else{0}
beta.n2[[m]][i,j]<- Judge
}
}}
result2 <- matrix(NA,100,3)
for (i in 1:3) {
for (j in 1:100) {
result2[j,i] <- mean(beta.n2[[i]][j,])
}
}
points(mu.X.true.seq,result2[,1],pch=3,lwd=1.5)
points(mu.X.true.seq,result2[,2],pch=3,col="red",lwd=1.5)
points(mu.X.true.seq,result2[,3],pch=3,col="blue",lwd=1.5)
Comments