[R] problem on simulation code (the loop unable to function effectively)

Jim Lemon drjimlemon at gmail.com
Tue Apr 19 02:25:32 CEST 2016


Hi Jeem,
First, please send questions like this to the help list, not me.

I assume that you are in a similar position to sjtan who has been
sending almost exactly the same questions.

The problem is not in the loops (which look rather familiar to me) but
in your initial assignments at the top. For instance:

scale parameter=(1,1.5,2,2.5,3)

produces an error which has nothing to do with the loops. This is a
very basic mistake, for:

scale_parameter<-c(1,1.5,2,2.5,3)

fixes it. I think if you learn a bit about basic R coding you will be
able to fix these problems yourself.

Jim


On Tue, Apr 19, 2016 at 4:03 AM,  <uk31429 at student.umt.edu.my> wrote:
> Greeting dr jim,
>
> I am student from Malaysia. I am doing R simulation study. For your
> information, I have been written a code relating to 2 gamma distribution
> with equal skewness.
> skewness=1.0
> shape parameter=16/9
> scale parameter=(1,1.5,2,2.5,3)
>
> Below are my coding, however, the code have some error and yet after
> trying this and that for a whole day, i couldn't spot the mistake. The
> output should be greater or smaller than 0.05 while no exceeding too much
> .But the for loop only able to function for scale parameter 1 .I try to
> apply another for loop for the scale parameter, but the output become
> worsen, as the simulation will even get hang. all of the output is 100.
> Please, could you give me some advice ?
>
> #For gamma disribution with equal skewness 1.5
> rm(list=ls())
> nSims<-100
> alpha<-0.05
>
> #here we declare the random seed generator
> set.seed(3)
>
> ## Put the samples sizes into matrix then use a loop for sample sizes
> sample_sizes<-matrix(c(10,10,10,25,25,25,25,50,25,100,50,25,50,100,100,25,100,100),nrow=2)
>
> #shape parameter for both gamma distribution for equal skewness
> shp<-rep(16/9,each=45)
>
> #scale parameter for sample 1
> #scale paramter for sample 2 set as constant 1
> d1<-matrix(c(1,1.5,2,2.5,3),ncol=1)
> scp<-rep(d1,9)
>
> #create a matrix combining the forty five cases of combination of sample
> sizes,shape and scale parameter
> all<- cbind(rep(sample_sizes[1,],5),rep(sample_sizes[2,],5),scp)
>
> # name the column samples 1 and 2 and standard deviation
> colnames(all) <- c("m", "n","scp")
>
>
> #set empty vector of length to store p-value
> equal3<-rep(0,nrow(all))
> unequal4<-rep(0,nrow(all))
> mann5<-rep(0,nrow(all))
>
> #set nrow =nsims because wan storing every p-value simulated
> #for gamma distribution with equal skewness
> matrix3_equal  <-matrix(0,nrow=nSims,ncol=3)
> matrix4_unequal<-matrix(0,nrow=nSims,ncol=3)
> matrix5_mann   <-matrix(0,nrow=nSims,ncol=3)
>
>
> # this loop steps through the all_combine matrix
>  for(ss in 1:nrow(all))
>
>   {  #generate samples from the first column and second column
>      m<-all[ss,1]
>      n<-all[ss,2]
>
>       for ( sim in 1:nSims)
>      {
>         #generate 2 random samples from gamma distribution with equal
> skewness
>         gamma1<-rgamma(m,16/9,all[ss,3])
>         gamma2<-rgamma(n,16/9,1)
>     #minus population mean from each sample to maintain the equality of
> null    #hypotheses (population mean =scale parameter *shape
> parameter)
> gamma1<-gamma1-16/9*all[ss,3]
> gamma2<-gamma2-16/9
>
>      matrix3_equal[sim,1]<-t.test(gamma1,gamma2,var.equal=TRUE)$p.value
>      matrix4_unequal[sim,2]<-t.test(gamma1,gamma2,var.equal=FALSE)$p.value
>      matrix5_mann[sim,3] <-wilcox.test(gamma1,gamma2)$p.value
>
>       }
> ##store the result
> equal3[ss]<- sum(matrix3_equal[,1]<alpha)
>    unequal4[ss]<-sum(matrix4_unequal[,2]<alpha)
>       mann5[ss]<- sum(matrix5_mann[,3]<alpha)
> }
> g<-cbind(all, equal3, unequal4, mann5)
>
> I will be really appreciated for any respond .Thanks.
>
> Your sincerely
> Jeem
>
>
>
>
> Computational Mathematics
> School of Informatics and Applied Mathematics
> Universiti Malaysia Terengganu
> 21030 Kuala Terengganu, Terengganu, Malaysia
> Phone : 017-7799039
> Email : uk31429 at student.umt.edu.my
>



More information about the R-help mailing list