Arrays and functions
On 13-Apr-08 08:43:33, Louisa Hay wrote:
Hi, I' am doing a stats project using R to work out the size of a t-test and wilcoxon test depending on the distribution and sample size. I just can't get it to work - I want to put my results from the function size() into an array. At the moment I keep getting the error message: Error in res[distribution, test, samplesize] <- results : subscript out of bounds Can anyone tell me where I'm going wrong,please?!
You have created a 3-dimensional array 'res' with
dimensions 2x2x3.
Apart from anything else, when you make the assignment
res[distribution,test,samplesize]<-results
the first value of 'samplesize' you use will be 10.
There is no array element with 3rd index equal to 10
(the maximum possible is 3). Hence (I suspect) your
error message.
I would recommend re-writing your code so that you use
the numbers of elements in each of your vectors
'distributions', 'tests' and 'samplesizes' to set the
ranges in your 'for' loops.
For example:
Ndist<-length(distributions)
Ntest<-length(tests)
Nsamp<-length(samplesizes)
and then you can do, for instance,
for(idist in (1:Ndist)){
for(itest in (1:Ntest)){
for(isamp in (1:Nsamp)){
for(i in (1:k)){
results<-size()
res[idist,itest,isamp]<-results
}
}
}
}
I also suspect that there are a lot of other things in
your code which could do with looking at too, but at
least the above approach should get you past the "subscript
out of bounds" problem!
Best wishes,
Ted.
size.power.test<-function(){
k<-1000
distributions<-c("Normal","Uniform")
tests<-c("t","Wilcoxon")
samplesizes<-c(10,30,40)
res<-array(0,c(length(distributions),
length(tests),
length(samplesizes))
)
dimnames(res)<-list(distributions,tests,samplesizes)
for(distribution in distributions){
for(test in tests){
for(samplesize in samplesizes){
for(i in 1:k){
results<-size()
res[distribution,test,samplesize]<-results
}
output<-size.power.test
return(output)
}
}
}
}
size<-function(k=1000,H0=0,mu.true=0,
sigma.true=1,alpha=0.05,
x1=-sqrt(3),y1=sqrt(3)){
distributions<-c("Normal","Uniform")
tests<-c("t","Wilcoxon")
samplesizes<-c(10,30,40)
reject1<-numeric(k)
for(distribution in distributions){
for(test in tests){
for(samplesize in samplesizes){
for(i in 1:k){
if(distribution=="Normal"){
#randomly generate data from normal distribution
#with true mean
dat<-rnorm(n=samplesize,mean=mu.true,sd=sigma.true)
}else{
#randomly generate data from a uniform distribution
dat<-runif(samplesize,min=x1,max=y1)
}
if(test=="t"){
#perform a t-test on sample
res<-t.test(dat,alternative="two.sided",mu=H0)
#record if p-value is less than alpha
reject1[i]<-(res$p.value<=alpha)
}else{
#perform a wilcoxon-test on sample
res<-wilcox.test(dat,alternative="two.sided",mu=H0)
#record if p-value is less than alpha
reject1[i]<-(res$p.value<=alpha)
}
}
output<-list(sum(reject1)/k)
return(output)
}
}
}
}
______________________________________________ R-help at 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.
-------------------------------------------------------------------- E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk> Fax-to-email: +44 (0)870 094 0861 Date: 13-Apr-08 Time: 10:23:11 ------------------------------ XFMail ------------------------------