在Python中随机生成特定长度的整数分区的算法?

klo*_*cey 6 python combinatorics sage

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

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

但是,我希望这对某人有帮助,我发现在函数返回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
Run Code Online (Sandbox Code Playgroud)

这增加了S找到长度分区的速率,并且看起来产生无偏的样本(我已经针对各种值的N和来检查了整个分区集的结果S).

然而,我正在使用N(例如10,000)和S(例如300)的值,这甚至使得这种方法不切实际地缓慢.与SAGE random_element()功能相关的评论承认有足够的优化空间.那么,有没有办法更快速地生成与给定值匹配的整数分区的无偏(即随机统一)样本,N并且S可能不生成不匹配的分区S?此外,在许多情况下使用共轭分区可以很好地生成无偏差的样本,但我不能说我完全理解为什么.

klo*_*cey 5

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

from sage.all import *
import random
Run Code Online (Sandbox Code Playgroud)

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

def min_max(n,s):

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

    return _min
Run Code Online (Sandbox Code Playgroud)

接下来,一个函数使用缓存和记忆化来找到 n 的分区数,其中 s 部分以 x 作为最大部分。这很快,但我认为有一个更优雅的解决方案。例如,通常:P(N,S,max=K) = P(NK,S-1) 感谢 ante ( /sf/users/34585351/ ) 帮助我解决这个问题: 找到数字给定总数、部分数量和最大被加数的整数分区

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)]
Run Code Online (Sandbox Code Playgroud)

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

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
Run Code Online (Sandbox Code Playgroud)