时间戳

首页 > 解决方案 > R中多个二项式随机数的模拟

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

问题描述

我有以下算法

Step 1. 生成 X1=x1~Bin(6,1/3)

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

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

调用函数 asmult_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,]

我没想到会有这样的结果。

为什么会有点？

我做错了什么？

标签： ralgorithmrandomstatisticssimulation

解决方案

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

这是我的解决方案。我的函数只在实验中返回。然后我使用复制返回多个（在我的情况下为 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 not 创建一个矩阵Matrix。

推荐阅读