Here is my code :
for (p in 1:10){
sample=rgamma(p,p,1)}
Here is what I get
[1] 1.841629
[1] 2.174076 1.410500
[1] 2.398601 4.674819 2.679830
[1] 2.736786 3.767747 4.546256 3.851677
[1] 4.204808 8.393887 10.312640 2.514957 4.863661
How can I simulate this for 1000 times and sum up the the results for each p.
Eg. for p=2, s=sum=(2.174076 1.410500 )
By this, there should be 1000 s value for each p
Thankyou very much.
CodePudding user response:
Using replicate
you can repeat a function for n
times.
simulate <- function(n) {
sample <- numeric(n)
for (p in 1:n){
sample[p] = sum(rgamma(p,p,1))
}
sample
}
res <- t(replicate(1000, simulate(10)))
CodePudding user response:
You can put your random gammas in an upper triangular matrix and sum the columns:
sz <- 1e3L
m <- matrix(0, nrow = sz, ncol = sz)
m[upper.tri(m, diag = TRUE)] <- rgamma((sz^2 sz)/2, rep(1:sz, 1:sz))
s <- colSums(m)
To verify that each column is getting the correct shape parameter:
sz <- 10L
m <- matrix(nrow = sz, ncol = sz)
m[upper.tri(m, diag = TRUE)] <- rep(1:sz, 1:sz)
> m
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 2 3 4 5 6 7 8 9 10
[2,] NA 2 3 4 5 6 7 8 9 10
[3,] NA NA 3 4 5 6 7 8 9 10
[4,] NA NA NA 4 5 6 7 8 9 10
[5,] NA NA NA NA 5 6 7 8 9 10
[6,] NA NA NA NA NA 6 7 8 9 10
[7,] NA NA NA NA NA NA 7 8 9 10
[8,] NA NA NA NA NA NA NA 8 9 10
[9,] NA NA NA NA NA NA NA NA 9 10
[10,] NA NA NA NA NA NA NA NA NA 10