将相同物品放入匿名桶后的可能分布及其概率

possible distributions and their probabilities after putting identical items into anonymous buckets

如果在其他地方很容易找到这个问题的答案,我们深表歉意。我的数学和统计数据很薄弱,因此我什至不知道我正在尝试做的事情的搜索词。 . .

我有 b 个匿名的不可区分的桶,我将 i 个相同的项目放入其中。我想知道所有可能的分布及其概率。例如,如果我有 3 个桶和 3 个项目,我想要的答案是:

请注意,存储桶是匿名的,因此我希望像上面那样合并相同的分布。比如[2,1,0]的情况其实就是[2,1,0],[0,2,1]等情况的总和

此外,我有最大存储桶大小的限制。比如3个球,3个桶,bucketsize=2应该return:

这可以通过概率树看出:

Insert item 1 into [0,0,0] -> [1,0,0] p=1  
Insert item 2 into [1,0,0] -> [2,0,0] p=1/3 OR [1,1,0] 2/3  
Insert item 3 into [2,0,0] -> [2,1,0] p=1.0   
Insert item 3 into [1,1,0] -> [2,1,0] p=2/3 OR [1,1,1] p=1/3  

So state [2,1,0] has two paths to it: 1/3*1 AND 2/3*2/3 = 7/9  
So state [1,1,1] has one path  to it: 2/3 * 1/3 = 2/9

这是每个概率的另一种观点: bins_and_balls http://bents.us/Pictures/bins_balls.png

非常感谢您阅读我的问题并提供任何帮助。我没有在 https://stats.stackexchange.com/ 上交叉发布这个,但如果人们认为它在那里更好,那么我会删除这个并在那里重新发布。

更新

评论中对一些建议算法的正确性进行了一些讨论。为了帮助验证,我编写了以下模拟器:

#! /usr/bin/env python

from __future__ import division
import random

def simulate(num_bucks,items,bsize,iterations=50000):
  perms = dict()
  for n in range(iterations):
    buckets = [0] * num_bucks 
    for i in range(items):
      while True:
        b = random.randint(0,num_bucks-1)
        if buckets[b] < bsize:
          break  # kludge, loop until we find an unfilled bucket 
      buckets[b] +=1
    buckets.sort()
    buckets = tuple(reversed(buckets))
    try:
      perms[buckets]['count'] += 1
    except KeyError:
      perms[buckets] = {'perm' : buckets, 'count' : 1}

  for perm in perms.values(): 
    perm['prob'] = perm['count'] / iterations
  return perms


def main():
  perms =    simulate(num_bucks=3,items=3,bsize=2)
  for perm in perms.values():
    print perm

if __name__ == "__main__":
    main()

其中有 3 个桶、3 个球、bucketsize 为 2 的输出,例如:

(1, 1, 1) 0.22394
(2, 1, 0) 0.77606

更新

我刚刚意识到您可能无法将我最初的建议用于较大的 n 和 r 值。至少,并非没有一些修改。例如,如果您将 4 件物品放入 5 个桶中,您可以(除其他外):

{3,1,0,0,0}
{2,2,0,0,0}

所以下面描述的 "n choose r" 技巧不起作用,因为有两种不同的方法可以只使用两个桶。

原回答

这是一个 "n choose r" 问题。组合。

对于 {3,0,0,0} 的情况,您问的是 "how many different ways can I choose one item from a set of four?" 答案是四。如果集合包含数字 1,2,3,4,那么您的组合是 {1},{2},{3},{4}.

对于 {2,1,0,0} 案例,您问的是从四个项目中选择两个项目有多少种不同的方法。答案是六个:{1,2},{1,3},{1,4},{2,3},{2,4},{3,4}.

而对于 {1,1,1,0} 的情况,答案是九。我就不一一列举了。

"n choose r"的计算结果是C(n,r) = n! / ( r! (n - r)! )

有关详细信息,请参阅 http://www.calculatorsoup.com/calculators/discretemathematics/combinations.php

概览

分发时,例如3 个项目超过 4 个桶,您可以将第一个项目放入 4 个桶中的任何一个,然后对于这 4 种可能性中的每一个,您可以将第二个项目放入 4 个桶中的任何一个,然后对于这 16 种可能性中的每一个,您可以将第三项放入 4 个桶中的任何一个,总共有 64 种可能性:

(3,0,0,0) (2,1,0,0) (2,0,1,0) (2,0,0,1) (2,1,0,0) (1,2,0,0) (1,1,1,0) (1,1,0,1) (2,0,1,0) (1,1,1,0) (1,0,2,0) (1,0,1,1) (2,0,0,1) (1,1,0,1) (1,0,1,1) (1,0,0,2) (2,1,0,0) (1,2,0,0) (1,1,1,0) (1,1,0,1) (1,2,0,0) (0,3,0,0) (0,2,1,0) (0,2,0,1) (1,1,1,0) (0,2,1,0) (0,1,2,0) (0,1,1,1) (1,1,0,1) (0,2,0,1) (0,1,1,1) (0,1,0,2) (2,0,1,0) (1,1,1,0) (1,0,2,0) (1,0,1,1) (1,1,1,0) (0,2,1,0) (0,1,2,0) (0,1,1,1) (1,0,2,0) (0,1,2,0) (0,0,3,0) (0,0,2,1) (1,0,1,1) (0,1,1,1) (0,0,2,1) (0,0,1,2) (2,0,0,1) (1,1,0,1) (1,0,1,1) (1,0,0,2) (1,1,0,1) (0,2,0,1) (0,1,1,1) (0,1,0,2) (1,0,1,1) (0,1,1,1) (0,0,2,1) (0,0,1,2) (1,0,0,2) (0,1,0,2) (0,0,1,2) (0,0,0,3)

您会注意到此列表中有很多重复项,例如(0,1,0,2) 出现 3 次。这是因为这 3 个项目是匿名且相同的,所以您可以通过多种方式将这 3 个项目分布到 4 个桶中并获得相同的结果,例如:

(0,0,0,0) → (0,1,0,0) → (0,1,0,1) → (0,1,0,2)
(0,0,0,0) → (0,0,0,1) → (0,1,0,1) → (0,1,0,2)
(0,0,0,0) → (0,0,0,1) → (0,0,0,2) → (0,1,0,2)

此外,由于存储桶是匿名且相同的,因此 (2,0,1,0) 和 (0,1,2,0) 等排列被认为是相等的。因此,这 64 种可能性组合在一起形成这 3 classes:

(3,0,0,0) x4
(2,1,0,0) x36
(1,1,1,0) x24

每个 class 的概率就是 class 的可能性的数量除以可能性的总数:

(3,0,0,0) → 4 / 64 = 1 / 16
(2,1,0,0) → 36 / 64 = 9 / 16
(1,1,1,0) → 24 / 64 = 6 / 16

生成所有可能性然后将它们组合到 classes 中的算法不是很有效,因为正如您所经历的那样,对于更多的项目和桶,可能性的数量很快就会变得巨大。
最好生成 classes(下面的第 1 步),计算每个 class 的排列数(下面的第 2 步),计算每个排列的重复分布数(下面的第 3 步),以及然后合并结果以获得每个 class 的概率(下面的第 4 步)。

第 1 步:生成分区

首先,生成 ipartitions (or more precisely: partitions with restricted number of parts) 和 b 部分:

(3,0,0,0) (2,1,0,0) (1,1,1,0)

生成分区的基本方法是使用递归算法,为第一部分选择所有可能的值,然后对其余部分进行递归,例如:

partition(4) :
• 4
• 3 + partition(1)
• 2 + partition(2)
• 1 + partition(3)
→ (4) (3,1) (2,2) (2,1,1) (1,3) (1,2,1) (1,1,2) (1,1,1,1)

为了避免重复,只生成非递增的值序列。为此,您必须告诉每个递归其部分的最大值是多少,例如:

partition(4, max=4) :
• 4
• 3 + partition(1, max=3)
• 2 + partition(2, max=2)
• 1 + partition(3, max=1)
→ (4) (3,1) (2,2) (2,1,1) (1,1,1,1)

限制部分数量的分区时,第一部分不应小于总值除以部分数量,例如将 4 分成 3 部分时,第一部分 1 没有分区,因为 1 < 4/3:

partition(4, parts=3, max=4) :
• 4,0,0
• 3 + partition(1, parts=2, max=3)
• 2 + partition(2, parts=2, max=2)
→ (4,0,0) (3,1,0) (2,2,0) (2,1,1)

当零件数为2时,这是最低的递归层次。然后通过迭代第一部分的值来生成分区,从最大值到值的一半,并将其余部分放在第二部分,例如:

partition(4, parts=2, max=3) :
• 3,1
• 2,2
→ (3,1) (2,2)

由于递归分区算法将每个任务减少为具有较小值和较少部分的多个任务,因此会多次生成具有较小值的相同分区。这就是为什么诸如记忆和制表之类的技术可以显着加快算法速度的原因。

第2步:计算排列数

然后,求每个分区的permutations个数:

(3,0,0,0) (0,3,0,0) (0,0,3,0) (0,0,0,3) → 4
(2,1,0,0) (2,0,1,0) (2,0,0,1) (1,2,0,0) (0,2,1,0) (0,2,0,1) (1,0,2,0) (0,1,2,0) (0,0,2,1) (1,0,0,2) (0,1,0,2) (0,0,1,2) → 12
(1,1,1,0) (1,1,0,1) (1,0,1,1) (0,1,1,1) → 4

您不需要生成所有这些排列来计算有多少;您只需检查步骤 1 中生成的每个分区中有多少重复值:

(3,0,0,0) (2,1,0,0) (1,1,1,0)

取桶数 b(4!)的阶乘,然后除以重复值数的阶乘:

(3,0,0,0) → 4! / (1! × 3!) = 4
(2,1,0,0) → 4! / (1! × 1! × 2!) = 12
(1,1,1,0) → 4! / (3! × 1!) = 4

第三步:计算分布数

然后,对于每个分区,我们必须找到可以在桶中分配项目的方式的数量,例如对于分区 (2,1,0,0):

• first item in bucket 1, second item in bucket 1, third item in bucket 2 → (2,1,0,0)
• first item in bucket 1, second item in bucket 2, third item in bucket 1 → (2,1,0,0)
• first item in bucket 2, second item in bucket 1, third item in bucket 1 → (2,1,0,0)

这个数字可以通过计算项目数量的阶乘 i (3!) 并将它除以值的阶乘的乘积(忽略零, 因为 0!= 1):

(3,0,0,0) → 3! / 3! = 1
(2,1,0,0) → 3! / (2! × 1!) = 3
(1,1,1,0) → 3! / (1! × 1! × 1!) = 6

第四步:计算概率

然后,对于每个分区,将排列数与分布数相乘,然后将它们相加得到可能性的总数:

4×1 + 12×3 + 4×6 = 64

或者简单地使用这个计算:

buckets ^ items = 4 ^ 3 = 64

然后每个分区的概率是排列数乘以分布数,除以可能性总数:

(3,0,0,0) → 4×1 / 64 = 4 / 64 = 1 / 16
(2,1,0,0) → 12×3 / 64 = 36 / 64 = 9 / 16
(3,0,0,0) → 4×6 / 64 = 24 / 64 = 6 / 16

代码示例

这是一个没有记忆的基本算法的例子:

function partitions(items, buckets) {
    var p = [];
    partition(items, buckets, items, []);
    probability();
    return p;

    function partition(val, parts, max, prev) {
        for (var i = Math.min(val, max); i >= val / parts; i--) {
            if (parts == 2) p.push({part: prev.concat([i, val - i])});
            else partition(val - i, parts - 1, i, prev.concat([i]));
        }
    }
    function probability() {
        var total = Math.pow(buckets, items);
        for (var i = 0; i < p.length; i++) {
            p[i].prob = permutations(p[i].part) * distributions(p[i].part) / total;
        }
    }
    function permutations(partition) {
        var dup = 1, perms = factorial(buckets);
        for (var i = 1; i < buckets; i++) {
            if (partition[i] != partition[i - 1]) {
                perms /= factorial(dup);
                dup = 1;
            } else ++dup;
        }
        return perms / factorial(dup);
    }            
    function distributions(partition) {
        var dist = factorial(items);
        for (var i = 0; i < buckets; i++) {
            dist /= factorial(partition[i]);
        }
        return dist;
    }
    function factorial(n) {
        var f = 1;
        while (n > 1) f *= n--;
        return f;
    }
}

var result = partitions(3,4);
for (var i in result) document.write(JSON.stringify(result[i]) + "<br>");
var result = partitions(7,5);
for (var i in result) document.write("<br>" + JSON.stringify(result[i]));


附加限制:存储桶大小

为了使算法给出的结果适应任何桶可以容纳的最大项目数的情况,您需要识别无效分布,然后在有效分布中重新分配它们的概率,例如:

i=5, b=4, s=5
5,0,0,0 → 1 / 256
4,1,0,0 → 15 / 256
3,2,0,0 → 30 / 256
3,1,1,0 → 60 / 256
2,2,1,0 → 90 / 256
2,1,1,1 → 60 / 256

i=5, b=4, s=4
4,1,0,0 → 16 / 256 ← (15 + 1)
3,2,0,0 → 30 / 256
3,1,1,0 → 60 / 256
2,2,1,0 → 90 / 256
2,1,1,1 → 60 / 256

i=5, b=4, s=3
3,2,0,0 → 35.333 / 256 ← (30 + 16 × 1/3)
3,1,1,0 → 70.666 / 256 ← (60 + 16 × 2/3)
2,2,1,0 → 90 / 256
2,1,1,1 → 60 / 256

i=5, b=4, s=2
2,2,1,0 → 172.444 / 256 ← (90 + 35.333 + 70.666 × 2/3)
2,1,1,1 → 83.555 / 256 ← (60 + 70.666 × 1/3)

下图显示了在将 7 项 3 桶分布限制为桶大小为 4 时必须重新计算哪些概率。

不受限制的存储桶大小的算法计算此图中的底部行,而无需创建整棵树。您可以调整它以直接计算任何单一分布的概率(请参阅下面的代码示例)。

所以你可以先运行原始算法,然后,对于每个无效分布,找出前一行中的哪些分布导致无效分布,并重新分配概率。

在图中所示的示例中,这意味着重新分配以红色表示的概率,这是 (4,0,0) 下的子树,其中 4 是最大桶大小。

所以一般来说,为了使原始算法找到的概率适应 s 的受限桶大小,您从 (s,0,0) 开始,找到此分布的所有子项:(s +1,0,0) 和 (s,1,0),计算它们的概率,取无效分布的概率并将它们重新分配到有效分布中,然后继续下一行,直到你重新计算了所有最后一行中具有值 s 的分布。

下面的代码示例,基于原始算法,展示了如何直接找到桶大小不受限制的任何分布的概率。这将是新算法中的一个基本步骤。

function probability(part) { // array of items per bucket
    var items = part.reduce(function(s, n) {return s + n}, 0);
    var buckets = part.length;
    var f = factorial(Math.max(items, buckets));
    var total = Math.pow(buckets, items);
    return permutations() * distributions() / total;

    function permutations() {
        var dup = 1, perms = f[buckets];
        for (var i = 1; i < buckets; i++) {
            if (part[i] != part[i - 1]) {
                perms /= f[dup];
                dup = 1;
            } else ++dup;
        }
        return perms / f[dup];
    }            
    function distributions() {
        return part.reduce(function(t, n) {return t / f[n]}, f[items]);
    }
    function factorial(max) {
        var f = [1];
        for (var i = 1; i <= max; i++) f[i] = f[i - 1] * i;
        return f;
    }
}

document.write(probability([2,2,1,1,1]));

以下代码示例对分区使用了不同的表示法:它使用的值是具有 index[ 的桶数,而不是值是每个桶的项目数的数组=177=] 项;所以 [0,3,2,0] 表示 3 个包含 1 个项目的桶和 2 个包含 2 个项目的桶,或常规表示法中的 [2,2,2,1,1]。

function probability(part) { // array of buckets per amount of items
    var items = part.reduce(function(t, n, i) {return t + n * i}, 0);
    var buckets = part.reduce(function(t, n) {return t + n}, 0);
    var f = factorial(Math.max(items, buckets));
    var total = Math.pow(buckets, items);
    return permutations() * distributions() / total;

    function permutations() {
        return part.reduce(function(t, n) {return t / f[n]}, f[buckets]);
    }            
    function distributions() {
        return part.reduce(function(t, n, i) {return t / Math.pow(f[i], n)}, f[items]);
    }
    function factorial(max) {
        var f = [1];
        for (var i = 1; i <= max; i++) f[i] = f[i - 1] * i;
        return f;
    }
}
document.write(probability([0,3,2,0,0,0,0,0]));

下面是使用替代符号的完整算法的实现;它生成有效分布(考虑桶大小),计算它们在不受限制的桶大小下的概率,然后重新分配无效分布的概率。

function Dictionary() {
    this.pair = {};
    this.add = function(array, value) {
        var key = array.join(',');
        if (key in this.pair) this.pair[key] += value;
        else this.pair[key] = value;
    }
    this.read = function() {
        for (var key in this.pair) {
            var prob = this.pair[key], array = key.split(',');
            array.forEach(function(v, i, a) {a[i] = parseInt(v);});
            delete this.pair[key];
            return {dist: array, prob: prob};
        }
    }
}
function partitions(items, buckets, size) {
    var dict = new Dictionary, init = [buckets], total = Math.pow(buckets, items), f = factorial();
    for (var i = 0; i < size; i++) init.push(0);
    partition(items, init, 0);
    if (size < items) redistribute();
    return dict;

    function partition(val, dist, min) {
        for (var i = min; i < size; i++) {
            if (dist[i]) {
                var d = dist.slice();
                --d[i]; ++d[i + 1];
                if (val - 1) partition(val - 1, d, i);
                else dict.add(d, probability(d));
            }
        }
        function probability(part) {
            return part.reduce(function(t, n, i) {return t / Math.pow(f[i], n);}, f[items])
            * part.reduce(function(t, n) {return t / f[n];}, f[buckets]) / total;
        }
    }
    function redistribute() {
        var parents = new Dictionary;
        --init[0]; ++init[size];
        parents.add(init, 0);
        for (var i = size; i < items; i++) {
            var children = new Dictionary;
            while (par = parents.read()) {
                var options = buckets - par.dist[size];
                par.prob += probability(par.dist, i) * par.dist[size] / buckets;
                for (var j = 0, b = 0; j < options; j += par.dist[b++]) {
                    while (par.dist[b] == 0) ++b;
                    var child = par.dist.slice();
                    --child[b]; ++child[b + 1];
                    if (i == items - 1) dict.add(child, par.prob * par.dist[b] / options);
                    else children.add(child, par.prob * par.dist[b] / options);
                }
            }
            parents = children;
        }
        function probability(part, items) {
            return part.reduce(function(t, n, i) {return t / Math.pow(f[i], n);}, f[items])
            * part.reduce(function(t, n) {return t / f[n];}, f[buckets]) / Math.pow(buckets, items);
        }
    }
    function factorial() {
        var f = [1], max = Math.max(items, buckets);
        for (var i = 1; i <= max; i++) f[i] = f[i - 1] * i;
        return f;
    }
}
var result = partitions(3,3,2);
for (var i in result.pair) document.write(i + " &rarr; " + result.pair[i] + "<br>");
var result = partitions(7,3,4);
for (var i in result.pair) document.write("<br>" + i + " &rarr; " + result.pair[i]);


(下面是第一个想法,后来发现效率太低了。)

在处理这个问题时,我想出了一个完全不同的算法版本,使用字典。它可能更接近您最初所做的,但应该更有效率。基本上,我使用了一种不同的方式来存储分布;我之前写的地方:

(3,2,2,1,0)

对于装有 3、2、2、1 和零项的 5 个桶,我现在使用:

(1,1,2,1)

对于 1 个空桶、1 个桶有 1 个项目、2 个桶有 2 个项目和 1 个桶有 3 个项目。此列表的长度由桶大小加 1 决定,不再由桶数决定。

这样做的好处是像这样的重复:

(3,2,2,1,0) (3,1,2,0,2) (1,2,2,3,0) (0,1,2,2,3) ...

现在所有符号都相同。

下面的代码示例使用了带有新符号的字典思想。 (不幸的是,事实证明它很慢,所以我会继续研究改进它。)

function Dictionary() {
    this.dict = [];
}
Dictionary.prototype.add = function(array, value) {
    var key = array.join(',');
    if (this.dict[key] == undefined) this.dict[key] = value;
    else this.dict[key] += value;
}
Dictionary.prototype.print = function() {
    for (var i in this.dict) document.write(i + " &rarr; " + this.dict[i] + "<br>");
}

function probability(items, buckets, size) {
    var dict = new Dictionary, dist = [buckets];
    for (var i = 0; i < size; i++) dist.push(0);
    p(dist, 0, 1);
    return dict;

    function p(dist, value, prob) {
        var options = buckets - dist[size];
        for (var i = 0, b = 0; i < options; i += dist[b++]) {
            while (dist[b] == 0) ++b;
            var d = dist.slice();
            --d[b]; ++d[b + 1];
            if (value + 1 < items) p(d, value + 1, prob * dist[b] / options);
            else dict.add(d, prob * dist[b] / options);
        }
    }
}

probability(3, 3, 3).print();
probability(3, 3, 2).print();
probability(5, 4, 5).print();
probability(5, 4, 2).print();