######################################################################
# R commands      Stat 704 
# Lab 4           University of South Carolina
# Simulation: Power
######################################################################
# This file contains the R commands for the lab.
#
# Lines beginning with the symbol '#' are comments in R.  All other
# lines contain code.
#
# In R for Windows, you may wish to open this file from the menu bar
# (File:Display file); you can then copy commands into the command
# window.  (Use the mouse to highlight one or more lines; then
# right-click and select "Paste to console".)
######################################################################



nosim<-10000
mu0<-30
mua<-32
sd<-4
n<-16
alt<-"greater"

SimData<-function(mua, sd, n){
	data<-rnorm(n=n, mean=mua, sd=sd)
	return(data)	
}

Calpvalue1<-function(data, mu0, alt){
	pvalue<-t.test(data, mu=mu0, alt=alt)$p.value
	return(pvalue)	
}
#######################################
# Version 1
#######################################

pvalues<-rep(NA, length=nosim)
for(i in 1:nosim){
	data<-SimData(mua=mua, sd=sd, n=n)
	pvalues[i]<-Calpvalue1(data, mu0=mu0, alt=alt)	
}
power<-mean(pvalues<0.05)
power

#######################################
# Version 2 
#######################################

Calpvalue2<-function(iter, mua, sd, n, mu0, alt){
	data<-SimData(mua=mua, sd=sd, n=n)
	pvalue<-t.test(data, mu=mu0, alt=alt)$p.value
	return(pvalue)	
}
pvalues<-sapply(1:nosim, Calpvalue2, mua=mua, sd=sd, n=n, mu0=mu0, alt=alt)

power<-mean(pvalues<0.05)
power