#pdf("C:\\normal_stuff.pdf")



set.seed(666)                        # set the seed for reproducible results
num.sample <- 10000                    # Number of samples
size.sample <- 100                   # Sample size for each sample
truemean <- 3                      # Assign the true mean of distribution to be sampled from
nullmean <- 3                      # Assign the mean under the null hypothesis
muhat <- rep(0,num.sample)        # Create vector to save mean of sample means
lrX2 <- rep(0,num.sample)         # Create vector to save LR Test Statistics
waldX2 <- rep(0,num.sample)       # Create vector to save Wald Test Statistics
scoreX2 <- rep(0,num.sample)      # Create vector to save Score Test Statistics

for (i in 1:num.sample) {
y <- rpois(size.sample,truemean)  # Generate size.sample Poisson(truemean) RVs
muhat[i] <- mean(y)               # Compute the mean of size.sample Poisson(truemean) RVs
lnLnull <- -size.sample*nullmean + sum(y)*log(nullmean)                      
# Compute log likelihood under null (ignoring terms not involving mu)
lnLmuhat <-  -size.sample*muhat[i] + sum(y)*log(muhat[i])                    
# Compute maximized log likelihood (ignoring terms not involving mu)
lrX2[i] <-  -2*(lnLnull - lnLmuhat)                      
# Compute Likelihood Ratio chi-square statistic
waldX2[i] <-  size.sample * ((muhat[i] - nullmean)^2 / muhat[i])                     
# Compute Wald chi-square statistic
score1  <- -size.sample + (sum(y) / nullmean)                          
# Compute Numerator of Score (LM) statistic
score2  <- size.sample / nullmean                         
# Compute Denominator of Score (LM) statistic
scoreX2[i] <- score1^2 / score2                    
# Compute Score chi-square statistic
}

(lr_power <- sum(lrX2 >= qchisq(0.95,1))/num.sample)           # Power of LR Test
(wald_power <- sum(waldX2 >= qchisq(0.95,1))/num.sample)       # Power of Wald Test
(score_power <- sum(scoreX2 >= qchisq(0.95,1))/num.sample)     # Power of Score Test

hist(muhat,main="",xlab="Mean Goals per Game",prob=T)
lines(density(muhat))



#dev.off()