【问题标题】:An algorithm for randomly generating integer partitions of a particular length, in Python?在 Python 中随机生成特定长度的整数分区的算法?
【发布时间】:2012-05-04 10:43:35
【问题描述】:

我一直在使用 SAGE 提供的 random_element() 函数为特定长度 (S) 的给定整数 (N) 生成随机整数分区。我正在尝试从给定值NS 的所有分区集中生成无偏随机样本。 SAGE 的函数快速返回 N 的随机分区(即Partitions(N).random_element())。

但是,添加S(即Partitions(N,length=S).random_element())时速度会大大降低。同样,过滤掉长度为SN 的随机分区非常慢。

但是,我希望这对某人有所帮助,我发现在函数返回 N 的分区与长度不匹配 S 的情况下,共轭分区的长度通常为 S。 :

S = 10
N = 100
part = list(Partitions(N).random_element())
    if len(part) != S:
        SAD = list(Partition(part).conjugate())
        if len(SAD) != S:
            continue

这增加了找到长度为S 的分区的速率,并且似乎产生了无偏样本(我已经针对NS 的各种值检查了整个分区集的结果)。

但是,我使用 N(例如 10,000)和 S(例如 300)的值,这使得即使这种方法也变得不切实际地慢。与 SAGE 的 random_element() 函数相关的评论承认有很大的优化空间。那么,有没有办法通过不生成与S 不匹配的分区来更快地生成匹配给定值NS 的整数分区的无偏(即随机统一)样本?此外,在许多情况下使用共轭分区可以很好地生成无偏样本,但我不能说我完全理解为什么。

【问题讨论】:

    标签: python combinatorics sage


    【解决方案1】:

    简单方法:随机分配整数:

    def random_partition(n, s):
        partition = [0] * s
        for x in range(n):
            partition[random.randrange(s)] += 1
        return partition
    

    【讨论】:

    • 感谢您的回复,但我看不到此函数如何根据均匀随机抽样产生分区。
    • @klocey,我错过了您从序列中生成随机元素的事实,抱歉。
    • 我实现了这个函数,并将它生成的随机样本与 N 和 S 的几种组合的完整分区集进行了比较。使用从分区方差生成的核密度曲线进行比较。与我尝试过的所有其他采样策略一样,此函数会产生有偏差的样本(低于预期方差的分区)。显然,对于给定的总 N 和长度 S,从所有分区的集合中生成无偏随机样本确实很困难。SAGE 函数是我所获得的最接近的函数,但它远非最优。
    【解决方案2】:

    我在尝试计算强生日问题的概率时遇到了类似的问题。

    首先,当只给出适量的数字时,分区函数会爆炸。你会返回很多信息。无论您使用哪种方法,N = 10000 和 S = 300 都会产生荒谬的数据量。它会很慢。您使用的任何纯 python 实现都有可能同样慢或慢。期待制作一个 CModule。

    如果您想尝试 python,我采用了结合 itertools 和生成器来降低内存使用率的方法。我的代码似乎不再方便了,但这是一个很好的实现:

    http://wordaligned.org/articles/partitioning-with-python

    编辑:

    找到我的代码:

    def partition(a, b=-1, limit=365):
      if (b == -1):
        b = a
      if (a == 2 or a == 3):
        if (b >= a and limit):
          yield [a]
        else:
          return
      elif (a > 3):
        if (a <= b):
          yield [a]
        c = 0
        if b > a-2:
          c = a-2
        else:
          c = b
        for i in xrange(c, 1, -1):
          if (limit):
            for j in partition(a-i, i, limit-1):
              yield [i] + j
    

    【讨论】:

    • 是的,组合爆炸很厉害。但是,我一次生成一个随机分区,并且只保留一个小的随机样本进行比较分析。我正在尝试为给定长度 S 的给定总 N 获取一个小的无偏随机分区样本。SAGE 的函数在 Cython 中运行,我自己的脚本也是如此,所以高效的速度并不像找到算法那么大或一种调整 SAGE 功能的方法,以避免生成不必要的分区(即长度不为 S 的分区)。我会看看你的实施和“强大的生日问题”。谢谢。
    • 找到了我的代码,它是一个生成器,它可以找到大小为 2 或更大的分区,最大为给定数量,您可以删除阻止分区小于 2 的逻辑。但我怀疑它会更快。
    【解决方案3】:

    最后,我有一个绝对无偏见的方法,拒绝率为零。当然,我已经对其进行了测试,以确保结果是整个可行集的代表性样本。它非常快而且完全没有偏见。享受吧。

    from sage.all import *
    import random
    

    首先,一个函数为 n 和 s 部分的分区找到最小的最大加数

    def min_max(n,s):
    
        _min = int(floor(float(n)/float(s)))
        if int(n%s) > 0:
            _min +=1
    
        return _min
    

    接下来,一个使用缓存和记忆来查找分区数量的函数 n 的 s 部分以 x 为最大部分。这很快,但我认为有 有一个更优雅的解决方案。例如,通常:P(N,S,max=K) = P(N-K,S-1) 感谢 ante (https://stackoverflow.com/users/494076/ante) 帮助我解决这个问题: Finding the number of integer partitions given a total, a number of parts, and a maximum summand

    D = {}
    def P(n,s,x):
        if n > s*x or x <= 0: return 0
        if n == s*x: return 1
        if (n,s,x) not in D:
            D[(n,s,x)] = sum(P(n-i*x, s-i, x-1) for i in xrange(s))
        return D[(n,s,x)]
    

    最后,一个函数找到 n 和 s 部分的均匀随机分区,没有拒绝率!每个随机选择的数字代码代表 n 个具有 s 个部分的特定分区。

    def random_partition(n,s):
        S = s
        partition = []
        _min = min_max(n,S)
        _max = n-S+1
    
        total = number_of_partitions(n,S)
        which = random.randrange(1,total+1) # random number
    
        while n:
            for k in range(_min,_max+1):
                count = P(n,S,k)
                if count >= which:
                    count = P(n,S,k-1)
                    break
    
            partition.append(k)
            n -= k
            if n == 0: break
            S -= 1
            which -= count
            _min = min_max(n,S)
            _max = k
    
        return partition
    

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2020-03-11
      • 2013-08-21
      • 2018-06-19
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多