Hi,I am currently doing a project in which we are to investigate the size and
power of three different one sample tests over three different distributions
using a number of different sample sizes and values for mu1. I have written a
function and was trying to get my answer for each test into the right position
in an array so the output is the power of each combination of test,
distribution, sample size and mu1, however my code is not working and i keep
getting an error message returned. My code is:
proj1.sim<-function(){
tests<-c("t","Wilcoxon","bootstrap")
distns<-c("Normal","Uniform","Beta")
sample.sizes<-c("5","10","25","50","100")
mu1s<-c("0","0.5","1","1.5","2")
res<-array(0,c(length(tests),length(distns),length(sample.sizes),length(mu1s)))
dimnames(res)<-list(tests,distns,sample.sizes,mu1s)
for(test in tests){
for(distn in distns){
for(sample.size in sample.sizes){
for(mu1 in mu1s){
result<-hypoth.test(sample.size,distn,test,mu1)
res[test,distn,sample.size,mu1]<-result
}}}}
output(res)
}
hypoth.test<-function(sample.size,distn,test,mu1){
#Purpose:
#Samples n values from the chosen distribution
#Inputs:
#n - sample size of data to generate
#distn - the distribution to generate data from
#test - which test to use
#H0 - the null hypothesis
#mu.true - the true value of the mean under the chosen distribution
#sigma.true - the true value of the variance under the chosen distribution
#alpha - trimming level
#K - number of simulations to perform
#Outputs:
#reject.null - the number of times we reject the null hypothesis
H0<-0
alpha<-0.05
K<-1100
sigma<-1
#create vectors for the results
reject<-numeric(K)
mu1<-numeric(5)
sample.size<-numeric(5)
#loop through the observations the required number of K simulations
for(i in 1:K){
#determine whether the data comes from a Normal or Uniform distribution
if(distn=="Normal"){
#Normal distribution
dat<-rnorm(n=sample.size,mean=mu1,sd=sigma)
} else if(distn=="Uniform"){
dat<-runif(n=sample.size,min=-sqrt(3)+mu1, max=sqrt(3)+mu1)
} else {
z<-rbeta(n=sample.size,shape1=0.8,shape2=7.2)
dat<-((z-0.1)*(sigma/0.1))+H0
}
#determine which test is required to sample the observations
if(test=="t"){
#t-test
res<-t.test(dat,alternative="two.sided",mu=H0)
} else if(test=="Wilcoxon"){
res<-wilcox.test(dat,alternative="two.sided",mu=H0)
} else {
nboot<-(K-1)
n<-length(dat)
bootmean<-NULL
for(i in 1:nboot){
rdata<-sample(dat,n,replace=T)
bootmean[i]<-mean(rdata)
}
nlo<-round((nboot+1)*(alpha/2))
nhi<-round((nboot+1)*((1-alpha)/2))
bootmean<-sort(bootmean)
low<-bootmean[nlo]
high<-bootmean[nhi]
res<-c(low,high)
}
#record whether we reject the null or not i.e. whether the p-value is less
#than or equal to the given alpha level(reject), or not (do not reject)
if(test=="bootstrap"){
reject[i]<-(H0<low|H0>high)
} else {
reject[i]<-(res$p.value<=alpha)
}
#sum together the total number of null hypothesis rejected
reject.null<-sum(reject)
}
power<-(reject.null)/K
return(power)
} Would you be able to tell me where I am going wrong? Thanks,Nicola Smith
_________________________________________________________________
[[elided Hotmail spam]]
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
