R：组合学，一种组合的排列数

Question

我有三个不同的事件 (1,2,3)，具有不同的概率 (0.15, 0.76, 0.09)，我想用替换画 5 次。 我现在可以使用

来确定可能的组合数量

nsimplex(3,5)  ### =21

来自组合包。 我可以使用

确定每个组合的概率

mySimplex <- xsimplex(3,5)
myProbs<-c(0.15, 0.76, 0.09)
results<- apply(mySimplex,2,dmultinom,prob=myProbs)

此外，我当然可以通过计算 3^5= 243 来确定排列的数量。

但是如果不手动计算它们，我怎么知道绘制相同组合的每个排列的频率？也就是说，我的每个组合中有多少个排列？

如果我没有正确理解，有 243 种排列构成 21 种不同的组合。现在我的问题是，每个组合有多少排列组合。例如。 {1,1,1,1,1} 组合只会被构建一次，而其他组合则由多个排列创建。

我想您可以通过使用每种组合的概率来实现这一点，但我不知道该怎么做？或者有没有其他方法可以在 R 中轻松确定？ 提前谢谢你。

Answer 1

第 1 项的 a 不可区分副本的排列数，第 2 项的 b，第 3 项的 c，其中 a + b + c = N 是 N! / (a! b! c!)。

例如，如果您有 (a,b,c) = (3,1,1)，则有 5!/(3! 1! 1!) = 20 安排。

  c  b  a  a  a    b  a  c  a  a    a  b  a  a  c    a  a  c  a  b    
  c  a  b  a  a    b  a  a  c  a    a  c  b  a  a    a  a  b  c  a    
  c  a  a  b  a    b  a  a  a  c    a  c  a  b  a    a  a  b  a  c    
  c  a  a  a  b    a  b  c  a  a    a  c  a  a  b    a  a  a  b  c    
  b  c  a  a  a    a  b  a  c  a    a  a  c  b  a    a  a  a  c  b  

一般情况下，我们可以这样计算数量

nperm<-function(...) {
   args<-as.numeric(list(...));
   num<-lfactorial(sum(args));
   den<-sum(lfactorial(args));
   return(round(exp(num-den)));
}

所以，例如，

x<-expand.grid(0:5,0:5,0:5)
x<-x[rowSums(x)==5,]
x[,"nperm"]<-apply(x,1,function(x) do.call(nperm,as.list(x)))

Var1 Var2 Var3 nperm
   5    0    0     1
   4    1    0     5
   3    2    0    10
   2    3    0    10
   1    4    0     5
   0    5    0     1
   4    0    1     5
   3    1    1    20
   2    2    1    30
   1    3    1    20
   0    4    1     5
   3    0    2    10
   2    1    2    30
   1    2    2    30
   0    3    2    10
   2    0    3    10
   1    1    3    20
   0    2    3    10
   1    0    4     5
   0    1    4     5
   0    0    5     1

正如预期的那样，还有sum(x[,"nperm"]) == 243。

Answer 2

为了使这个可重现，我需要使用set.seed(<some_value>)，但这是使用sample 绘制不同组合的一次尝试（不考虑排列不同。如果要考虑排列不同，则取出sort 步骤：

table(                    # get the counts of distinct combinations
   apply(                 # this will collapse values by column
      replicate(100000,   # yields a 100,000 column matrix
             {sample(c("1","2","3"), 5 ,repl=TRUE, prob=c(.5,.25,.25) )}), 
           2, function(x) paste(sort(x), collapse=".")) )

1.1.1.1.1 1.1.1.1.2 1.1.1.1.3 1.1.1.2.2 1.1.1.2.3 1.1.1.3.3 1.1.2.2.2 
     3090      7705      8144      7851     15408      7649      3997 
1.1.2.2.3 1.1.2.3.3 1.1.3.3.3 1.2.2.2.2 1.2.2.2.3 1.2.2.3.3 1.2.3.3.3 
    11731     11554      3940       949      3844      5955      4019 
1.3.3.3.3 2.2.2.2.2 2.2.2.2.3 2.2.2.3.3 2.2.3.3.3 2.3.3.3.3 3.3.3.3.3 
      961        99       506       990       997       510       101 

A.Webb 建议我们将理论 dmultinom 与实践进行比较：

dmultinom(c(4,1,0),prob=c(0.5,0.25,0.25))*2
[1] 0.15625

因此，第一个值 3125 的预测看起来可以说是准确的，而不是在 3090 处模拟的，第二个和第三个值也是 7812.5 与 7705 和 8144 的预测。