0

i have combined all the variables in a matrix, and i wish to conduct a simulation row by row. But i found out the code only works for nine samples only ,and not all of them (45). I tried , the loop is iterating but because of these two lines , so problem occurs.

 #minus the population mean to ensure the true of null hypo 
            gamma1<-gamma1-16/9*all_combine1[ss,4]
            gamma2<-gamma2-16/9

can anyone pls help ....please

 #For gamma disribution with equal skewness 1.5

#to evaluate the same R function on many different sets of data
library(parallel)

nSims<-100   #number of simulation
alpha<-0.05   #significance level

#set nrow =nsims because wan storing every p-value simulated
matrix3_equal  <-matrix(0,nrow=nSims,ncol=3)
matrix4_unequal<-matrix(0,nrow=nSims,ncol=3)
matrix5_mann   <-matrix(0,nrow=nSims,ncol=3)
#set empty vector of length to store p-value 
equal3<-c(rep(0,nrow(all_combine1)))
unequal4<-c(rep(0,nrow(all_combine1)))
mann5<-c(rep(0,nrow(all_combine1)))

#for gamma distribution with equal skewness
# to ensure the reproducity of the result 
#here we declare the random seed generator
set.seed(1)

## 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
scp1<-c(1,1.5,2,2.5,3)
scp1<-rep(scp1,9)   

#create a matrix combining the forty five cases of combination of sample sizes,shape and scale parameter
all_combine1 <- cbind(rep(sample_sizes[1,], 5),rep(sample_sizes[2,],5),shp,scp1)

# name the column samples 1 and 2 and standard deviation
colnames(all_combine1) <- c("m", "n","sp(skewness1.5)","scp1")

##for the samples sizes into matrix then use a loop for sample sizes
 # this loop steps through the all_combine matrix
  for(ss in 1:nrow(all_combine1))  
  {
    #generate samples from the first column and second column
     m<-all_combine1[ss,1]
     n<-all_combine1[ss,2]   

       for (sim in 1:nSims)
       {
        #generate 2 random samples from gamma distribution with equal skewness
        gamma1<-rgamma(m,16/9,all_combine1[ss,4])
        gamma2<-rgamma(n,16/9,1)

      #minus the population mean to ensure the true of null hypo 
        gamma1<-gamma1-16/9*all_combine1[ss,4]
        gamma2<-gamma2-16/9

        #extract p-value out and store every p-value into matrix
        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)

}

This is my result.Clearly it do not run successfully for every loop.TT

        m   n sp(skewness1.5) scp1 equal3 unequal4 mann5
 [1,]  10  10        1.777778  1.0      9        9     6
 [2,]  10  25        1.777778  1.5     94       93    95
 [3,]  25  25        1.777778  2.0    100      100   100
 [4,]  25  50        1.777778  2.5    100      100   100
 [5,]  25 100        1.777778  3.0    100      100   100
 [6,]  50  25        1.777778  1.0      3        8     6
 [7,]  50 100        1.777778  1.5    100      100   100
 [8,] 100  25        1.777778  2.0    100      100   100
 [9,] 100 100        1.777778  2.5    100      100   100
[10,]  10  10        1.777778  3.0    100      100   100
[11,]  10  25        1.777778  1.0      3        4     6
[12,]  25  25        1.777778  1.5     99       99   100
[13,]  25  50        1.777778  2.0    100      100   100
[14,]  25 100        1.777778  2.5    100      100   100
[15,]  50  25        1.777778  3.0    100      100   100
[16,]  50 100        1.777778  1.0      3        4     1
[17,] 100  25        1.777778  1.5    100      100   100
[18,] 100 100        1.777778  2.0    100      100   100
[19,]  10  10        1.777778  2.5    100      100   100
[20,]  10  25        1.777778  3.0    100      100   100
[21,]  25  25        1.777778  1.0      4        3     5
[22,]  25  50        1.777778  1.5    100       99   100
[23,]  25 100        1.777778  2.0    100      100   100
[24,]  50  25        1.777778  2.5    100      100   100
[25,]  50 100        1.777778  3.0    100      100   100
[26,] 100  25        1.777778  1.0      8        9    10
[27,] 100 100        1.777778  1.5    100      100   100
[28,]  10  10        1.777778  2.0    100      100   100
[29,]  10  25        1.777778  2.5    100      100   100
[30,]  25  25        1.777778  3.0    100      100   100
[31,]  25  50        1.777778  1.0      2        3     2
[32,]  25 100        1.777778  1.5    100      100   100
[33,]  50  25        1.777778  2.0    100      100   100
[34,]  50 100        1.777778  2.5    100      100   100
[35,] 100  25        1.777778  3.0    100      100   100
[36,] 100 100        1.777778  1.0      7        7     5
[37,]  10  10        1.777778  1.5     88       87    90
[38,]  10  25        1.777778  2.0    100      100   100
[39,]  25  25        1.777778  2.5    100      100   100
[40,]  25  50        1.777778  3.0    100      100   100
[41,]  25 100        1.777778  1.0      7        7     6
[42,]  50  25        1.777778  1.5    100      100   100
[43,]  50 100        1.777778  2.0    100      100   100
[44,] 100  25        1.777778  2.5    100      100   100
[45,] 100 100        1.777778  3.0    100      100   100
j.l
  • 25
  • 8
  • Just noticed that only scale parameter 1 is being running!!! Trying to figure out why – j.l Apr 18 '16 at 10:27

2 Answers2

1

I believe your error lies with these lines:

  ##store the result
  equal[ss]<- mean(matrix2_equal[,1]<=alpha)
  unequal[ss]<-mean(matrix5_unequal[,2]<=alpha)
  mann[ss]<- mean(matrix8_mann[,3]<=alpha)

matrix2_equal[,1]<=alpha will return a value of true or false thus mean(matrix2_equal[,1]<=alpha) is basically returning the % of the True from your model. This might be what you wanted:

equal[ss]<- mean(matrix2_equal[matrix2_equal[,1]<=alpha, 1])

FYI: This question is related to your previous post of: R： coding why show 0.00 in result

Community
  • 1
  • 1
Dave2e
  • 22,192
  • 18
  • 42
  • 50
  • i have to admit that there are related to my previous question . Sorry for that ....why in this case , it will unable to function ? before this, i try to apply same idea for simulation on normal distribution, it works without adding the part that you suggested... I tried yours suggestion , it works , i get the result 1.795561e-02 (example of one of it) like this. – j.l Apr 18 '16 at 02:56
  • I would have look the simulation some more, but I am guessing because the normal distribution is symmetric the function provide the correct answer with the incorrect calculation and since gamma is skewed, it showed a bias. If this answer worked, please accept the answer to close the question. - Thanks. – Dave2e Apr 18 '16 at 03:03
  • can i ask further about the purpose of adding the [matrix2_equal[,1]<=alpha, 1] into the coding...? isnt it mean that we are taking the element obtain from matrix2_equal[,1]<=alpha ? – j.l Apr 18 '16 at 03:07
  • it do works , the answer do changes but the answer is not correct as what it should be for example, it cannot be this numbers 215147e-51 for type1error ...i am sorry but i will think for another way >< 2 – j.l Apr 18 '16 at 03:25
  • What is the purpose of: gamma1<-gamma1-16/9*all_combine1[ss,4] it appears you are over correcting gamma1 with respect to gamma2. I suggest commenting out these lines and see if that improves your results. – Dave2e Apr 18 '16 at 14:35
  • These two lines is to ensure the equality of population mean ,because standard deviation is dependent of mean, thus in order to have equality of mean, have to minus the mean from each sample. i tried to taking out the two lines ,but it only help to getting few results out... – j.l Apr 18 '16 at 15:10
1

@Dave2e commenting out the two lines give the following result.

 m   n scp equal3 unequal4 mann5
 [1,]  10  10 1.0      8        8     9
 [2,]  10  25 1.5     16       36    23
 [3,]  25  25 2.0     83       82    78
 [4,]  25  50 2.5    100      100   100
 [5,]  25 100 3.0    100      100   100
 [6,]  50  25 1.0      3        5     7
 [7,]  50 100 1.5     82       86    79
 [8,] 100  25 2.0     98       92    91
 [9,] 100 100 2.5    100      100   100
[10,]  10  10 3.0     76       72    77
[11,]  10  25 1.0      1        3     3
[12,]  25  25 1.5     44       42    37
[13,]  25  50 2.0     94       96    92
[14,]  25 100 2.5    100      100   100
[15,]  50  25 3.0    100      100   100
[16,]  50 100 1.0      4        4     3
[17,] 100  25 1.5     72       54    56
[18,] 100 100 2.0    100      100   100
[19,]  10  10 2.5     65       60    57
[20,]  10  25 3.0     90       98    95
[21,]  25  25 1.0      2        2     5
[22,]  25  50 1.5     48       61    50
[23,]  25 100 2.0     95       96    93
[24,]  50  25 2.5    100       99    98
[25,]  50 100 3.0    100      100   100
[26,] 100  25 1.0      5        6     2
[27,] 100 100 1.5    100      100    95
[28,]  10  10 2.0     50       49    49
[29,]  10  25 2.5     79       92    85
[30,]  25  25 3.0     99       99    99
[31,]  25  50 1.0      6        3     6
[32,]  25 100 1.5     58       76    54
[33,]  50  25 2.0     94       91    90
[34,]  50 100 2.5    100      100   100
[35,] 100  25 3.0    100      100   100
[36,] 100 100 1.0      3        3     4
[37,]  10  10 1.5     22       20    13
[38,]  10  25 2.0     45       70    55
[39,]  25  25 2.5     97       97    95
[40,]  25  50 3.0    100      100   100
[41,]  25 100 1.0      5        5     5
[42,]  50  25 1.5     62       48    52
[43,]  50 100 2.0    100      100   100
[44,] 100  25 2.5    100      100   100
[45,] 100 100 3.0    100      100   100
j.l
  • 25
  • 8