Monte Carlo Simulation for Power Function

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 a 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, and only the values which are abnormally large will be rejected.

With such a setting, if we let D = 1 when rejecting H0, and D = 0 if not rejected, 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 the critical value, we are trying to find the tail probability of a series of normal distributions given the 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 the 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 an estimated sigma does not have such a good fit.

R Code (click to watch)

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)

Page updated

Google Sites

Report abuse