r - R中多个二项式随机数的模拟

标签 r algorithm random statistics simulation

我有以下算法

第一步,生成X1=x1~Bin(6,1/3)

步骤 2. 生成 X2|X1=x1~Bin(6-x1,(1/3)/(1-1/3))

第三步生成X3|X1=x1,X2=x2~Bin(6-x1-x2,(1/3)/(1-1/3-1/3))

第4步。重复第1-3步N次。

这是我在 R 中实现此算法的方法:

mult_binom<-function(n) #n=6
{
  n=1000
  random_vectors<-Matrix(0,n,3)

  for(i in 1:n){

    X1<-rbinom(n,3,1/3) 

    X2<-rbinom(n-X1,3,(1/3)/(1-(1/3)))

    X3<-rbinom(n-X1-X2,3,(1/3)/(1-(1/3)-(1-3)))


    arr<-c(X1,X2,X3)


  }
  for(j in 1:n){

    random_vectors[j]<-arr[j]
  }
  return(random_vectors)
}

将函数调用为 mult_bin(6) 会产生一个类似的矩阵,如下所示

1000 x 3 sparse Matrix of class "dgCMatrix"

  [1,] 1 . .
  [2,] 1 . .
  [3,] 1 . .
  [4,] 2 . .
  [5,] 1 . .
  [6,] 1 . .
  [7,] 1 . .
  [8,] . 3 .

一直到 [1000,]

没想到会是这样的结果。

为什么会有点?

我做错了什么?

最佳答案

您的实现中有几个错误。最重要的是,rbinom 的第一个参数不是二项分布的参数n,而是您要生成的随机数的个数。

这是我的解决方案。我的函数只在实验中返回。然后我使用 replicate 返回多个(在我的例子中是 5 个)实验的结果:

myfun <- function(){

  x1 <- rbinom(1, 6, 1/3) 
  x2 <- rbinom(1, 6 - x1, (1/3)/(1-(1/3)))
  x3 <- rbinom(1, 6 - x1 - x2, (1/3)/(1-(1/3)-(1/3)))

  return(c(X1 = x1, X2 = x2, X3 = x3))  
}

    set.seed(1)
    replicate(5, myfun())
   [,1] [,2] [,3] [,4] [,5]
X1    1    4    4    0    3
X2    2    0    1    2    1
X3    3    2    1    4    2

在此输出中,每一列都是一个实验的结果。你可以看到数字总和为 6。 另请注意,我使用 set.seed 设置了一个随机种子。这可确保您的结果可重现。

在您的输出中出现这些点是因为您使用 Matrix 包创建了一个 Matrix 对象,而不是使用“普通”矩阵。通常,您使用 matrix 而不是 Matrix 创建矩阵。

关于r - R中多个二项式随机数的模拟,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/55763855/

相关文章:

r - 合并具有多个匹配项的数据框时仅选择第一行

在 R 中按组删除异常值

r - 取消列出列以在数据框中创建唯一行

r - 使用 dplyr 和 group_by 编写自己的函数 - 如何继续更改列名

用于从每个元素具有不同概率的列表中进行选择的 C++ 函数

c++ - C++ 中简单的类型驱动随机模型构造

c - 播种许多 PRNG,然后不得不再次播种,什么是高质量的方法?

c# - Merge 2排序时间序列算法

arrays - 总和等于键的数组的最小子集。条件 : Values can be used any number of times

algorithm - 将 M 人分成 N 个团队,并设置比例限制