r - 如何多次模拟不同的样本量并计算每个样本量的总和

Question

这是我的代码：

for (p in 1:10){
sample=rgamma(p,p,1)}

这是我得到的

[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

如何模拟 1000 次并总结每个 p 的结果。
例如。对于 p=2, s=sum=(2.174076 + 1.410500 )
这样，每个 p 应该有 1000 s 值
非常感谢。

Answer 1

你可以把你的随机伽马放在一个上三角矩阵中并对列求和：

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)

要验证每列是否获得正确的形状参数：

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

编辑：或稍快的单线：

s <- colSums(replace(matrix(0, sz, sz), sequence(1:sz, seq(1, sz^2 - sz + 1, sz)), rgamma((sz^2 + sz)/2, rep(1:sz, 1:sz))))

或者您可以直接对总和进行抽样：

s <- rgamma(sz, (1:sz)^2)

Answer 2

使用replicate您可以n多次重复一个功能。

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)))