2

我们如何重新编码一组严格递增(或严格递减)的正整数P,以减少整数之间可能出现的正整数的数量?将一组大整数变换成一组小整数

我们为什么要这样做:假设我们想随机抽样P但是1.)P太大而无法枚举,并且2.)P的成员以非随机的方式相关,但在某种程度上这太复杂了,无法进行采样。但是,当我们看到它时,我们就知道P的一个成员。假设我们知道P [0]和P [n],但是不能理解枚举全部P或者理解P的成员是如何相关的。同样,在P [0]和P [n]之间出现的所有可能整数的数量比P的大小大很多倍,这使得随机抽取P成员的机会非常小。

例:让P [0] = 2101010101 & P [N] = 505050505.现在,也许我们只在P [0],P [N]具有特定质量之间的整数兴趣(如P [x]中的所有整数小于或等于Q,P的每个成员的最大整数为7或更小)。因此,并非所有正整数P [n] < = X < = P [0]属于P.我感兴趣的P在下面的评论中讨论。

我已经试过什么:如果P是一个严格递减的一套,我们知道P [0],P [N],那么我们可以把每一个成员,就好像它是与P减去[0]。这样做会减少每个数字,也许会大大减少并将每个成员保持为唯一的整数。对于我感兴趣的P(下),可以将P的每个减少的值除以一个公分母(9,11,99),这减少了P成员之间可能的整数数目。发现结合使用时,这些方法将所有P [0]的集合减少了几个数量级,使得随机从所有正整数P [n]中抽取一个P成员的机会] < = X < = P [0]仍然非常小。

注意:应该清楚,我们必须知道关于体育的一些事情。如果我们不这样做,那基本上意味着我们不知道我们在找什么。当我们在P [0]和P [n]之间随机抽样整数(记录或不记录)时,我们需要能够说“Yup,属于P.”,如果确实如此。

好的答案可以大大增加我开发的计算算法的实际应用。我对此感兴趣的一种P的例子在评论2中给出。我坚持要给予应有的信贷。

+0

你有能力将k映射到P [k]吗?你怎么知道P [k]是什么?其中一些可能取决于这些P [k]实际是什么以及它们是如何表示的。 – mhum 2013-02-16 00:02:40

+0

P的集合中的每个成员对应于Q的整数分区,N部分。例如令Q = 20且N = 4,则对于Q和N的第一词汇划分,即[17,1,1,1]映射到P,为17010101(即P [0])。因此,整数分区[5,5,5,5]对应于5050505。从具有N个部分的Q的所有整数分区集合中随机抽样,无递归且具有低的拒绝率。在那里,我把它全部拿走了。为了答案,这是值得的。 – klocey 2013-02-16 00:11:48

+0

不妨加上一些鼓励:当N小于5时,这种方法快速致盲。快速致盲,我的意思是它使事件的概率太小,以至于我的计算机无法计算它(例如,使用典型的随机分区算法)变得非常可能。 – klocey 2013-02-16 00:14:42

回答

1

尽管最初的问题是询问关于整数编码的非常通用的场景,但我会建议不太可能存在一种完全通用的方法。例如,如果P [i]或多或少是随机的(从信息论的角度来看),如果有什么工作的话,我会感到惊讶。

因此,让我们把注意力转向OP生成包含完全K部分的整数N分区的实际问题。当用组合对象作为整数进行编码时,我们应该尽可能多地保留组合结构。 为此,我们转向经典文本Combinatorial Algorithms by Nijenhuis and Wilf,特别是第13章。实际上,在本章中,他们演示了一个框架来枚举和抽样一些组合族 - 包括N的分区,其中最大的部分等于K.使用众所周知的分区与K部分之间的分区和最大部分为K的分区(采用Ferrers diagram的转置),我们发现我们只需要对解码过程进行更改。

不管怎么说,这里的一些源代码:

import sys 
import random 
import time 

if len(sys.argv) < 4 : 
    sys.stderr.write("Usage: {0} N K iter\n".format(sys.argv[0])) 
    sys.stderr.write("\tN = number to be partitioned\n") 
    sys.stderr.write("\tK = number of parts\n") 
    sys.stderr.write("\titer = number of iterations (if iter=0, enumerate all partitions)\n") 
    quit() 

N = int(sys.argv[1]) 
K = int(sys.argv[2]) 
iters = int(sys.argv[3]) 

if (N < K) : 
    sys.stderr.write("Error: N<K ({0}<{1})\n".format(N,K)) 
    quit() 

# B[n][k] = number of partitions of n with largest part equal to k 
B = [[0 for j in range(K+1)] for i in range(N+1)] 

def calc_B(n,k) : 
    for j in xrange(1,k+1) : 
     for m in xrange(j, n+1) : 
      if j == 1 : 
       B[m][j] = 1 
      elif m - j > 0 : 
       B[m][j] = B[m-1][j-1] + B[m-j][j] 
      else : 
       B[m][j] = B[m-1][j-1] 

def generate(n,k,r=None) : 
    path = [] 
    append = path.append 

    # Invalid input 
    if n < k or n == 0 or k == 0: 
     return [] 

    # Pick random number between 1 and B[n][k] if r is not specified 
    if r == None : 
     r = random.randrange(1,B[n][k]+1) 

    # Construct path from r  
    while r > 0 : 
     if n==1 and k== 1: 
      append('N') 
      r = 0 ### Finish loop 
     elif r <= B[n-k][k] and B[n-k][k] > 0 : # East/West Move 
      append('E') 
      n = n-k 
     else : # Northeast/Southwest move 
      append('N') 
      r -= B[n-k][k] 
      n = n-1 
      k = k-1 

    # Decode path into partition  
    partition = [] 
    l = 0 
    d = 0  
    append = partition.append  
    for i in reversed(path) : 
     if i == 'N' : 
      if d > 0 : # apply East moves all at once 
       for j in xrange(l) : 
        partition[j] += d 
      d = 0 # reset East moves 
      append(1) # apply North move 
      l += 1    
     else : 
      d += 1 # accumulate East moves  
    if d > 0 : # apply any remaining East moves 
     for j in xrange(l) : 
      partition[j] += d 

    return partition 


t = time.clock() 
sys.stderr.write("Generating B table... ")  
calc_B(N, K) 
sys.stderr.write("Done ({0} seconds)\n".format(time.clock()-t)) 

bmax = B[N][K] 
Bits = 0 
sys.stderr.write("B[{0}][{1}]: {2}\t".format(N,K,bmax)) 
while bmax > 1 : 
    bmax //= 2 
    Bits += 1 
sys.stderr.write("Bits: {0}\n".format(Bits)) 

if iters == 0 : # enumerate all partitions 
    for i in xrange(1,B[N][K]+1) : 
     print i,"\t",generate(N,K,i) 

else : # generate random partitions 
    t=time.clock() 
    for i in xrange(1,iters+1) : 
     Q = generate(N,K) 
     print Q 
     if i%1000==0 : 
      sys.stderr.write("{0} written ({1:.3f} seconds)\r".format(i,time.clock()-t)) 

    sys.stderr.write("{0} written ({1:.3f} seconds total) ({2:.3f} iterations per second)\n".format(i, time.clock()-t, float(i)/(time.clock()-t) if time.clock()-t else 0)) 

而这里的性能的一些例子(在MacBook Pro上8.3,2GHz的酷睿i7,4 GB,Mac OSX版10.6.3,Python的2.6.1):

mhum$ python part.py 20 5 10 
Generating B table... Done (6.7e-05 seconds) 
B[20][5]: 84 Bits: 6 
[7, 6, 5, 1, 1] 
[6, 6, 5, 2, 1] 
[5, 5, 4, 3, 3] 
[7, 4, 3, 3, 3] 
[7, 5, 5, 2, 1] 
[8, 6, 4, 1, 1] 
[5, 4, 4, 4, 3] 
[6, 5, 4, 3, 2] 
[8, 6, 4, 1, 1] 
[10, 4, 2, 2, 2] 
10 written (0.000 seconds total) (37174.721 iterations per second) 

mhum$ python part.py 20 5 1000000 > /dev/null 
Generating B table... Done (5.9e-05 seconds) 
B[20][5]: 84 Bits: 6 
100000 written (2.013 seconds total) (49665.478 iterations per second) 

mhum$ python part.py 200 25 100000 > /dev/null 
Generating B table... Done (0.002296 seconds) 
B[200][25]: 147151784574 Bits: 37 
100000 written (8.342 seconds total) (11987.843 iterations per second) 

mhum$ python part.py 3000 200 100000 > /dev/null 
Generating B table... Done (0.313318 seconds) 
B[3000][200]: 3297770929953648704695235165404132029244952980206369173 Bits: 181 
100000 written (59.448 seconds total) (1682.135 iterations per second) 

mhum$ python part.py 5000 2000 100000 > /dev/null 
Generating B table... Done (4.829086 seconds) 
B[5000][2000]: 496025142797537184410324290349759736884515893324969819660 Bits: 188 
100000 written (255.328 seconds total) (391.653 iterations per second) 

mhum$ python part-final2.py 20 3 0 
Generating B table... Done (0.0 seconds) 
B[20][3]: 33 Bits: 5 
1 [7, 7, 6] 
2 [8, 6, 6] 
3 [8, 7, 5] 
4 [9, 6, 5] 
5 [10, 5, 5] 
6 [8, 8, 4] 
7 [9, 7, 4] 
8 [10, 6, 4] 
9 [11, 5, 4] 
10 [12, 4, 4] 
11 [9, 8, 3] 
12 [10, 7, 3] 
13 [11, 6, 3] 
14 [12, 5, 3] 
15 [13, 4, 3] 
16 [14, 3, 3] 
17 [9, 9, 2] 
18 [10, 8, 2] 
19 [11, 7, 2] 
20 [12, 6, 2] 
21 [13, 5, 2] 
22 [14, 4, 2] 
23 [15, 3, 2] 
24 [16, 2, 2] 
25 [10, 9, 1] 
26 [11, 8, 1] 
27 [12, 7, 1] 
28 [13, 6, 1] 
29 [14, 5, 1] 
30 [15, 4, 1] 
31 [16, 3, 1] 
32 [17, 2, 1] 
33 [18, 1, 1] 

我将它留给OP来验证此代码确实根据所需的(统一)分布生成分区。

编辑:增加了一个枚举功​​能的例子。

+0

我将函数的输出与已知无偏但较慢的函数的输出进行了比较。我生成了每个函数的随机样本,并比较了goo.gl/uNWR4中统计特征的分布(偏度,均匀度,中位数和值)。如果函数没有偏差,则所得到的核密度曲线不会显示任何系统差异。事实上,它们很接近但不同。所以,一种偏见已经蔓延开来。我将在明天花费你研究你所做的事情(这很酷),并试图找出偏见在哪里。如果你喜欢,我可以通过电子邮件发布/输出结果/脚本。 – klocey 2013-02-18 03:11:43

+0

好的。我只测试了一些N&K足够小的情况,以便我可以为所有可能的分区获得体面的样本。你能描述这些差异吗?你正在测试哪个N&K?我会再试一次,并尝试调试。 – mhum 2013-02-18 03:18:35

+0

我使用N = {20,40}&K = {3,5,10}足够小来检查可行集合(即N&K的所有分区)的分布特征,我使用核密度曲线比较了随机样本生成的w /你的方法。统计均匀度高于应有的水平。偏度往往低于应有的水平。中等加数值显示多模态,但应该是单峰的。随机样本的K密度曲线彼此高度一致,但与可行集不一致。 – klocey 2013-02-18 09:07:57

0

下面是一个脚本,它完成了我所问的内容,至于重新编码整数表示用K个部分表示N的整数分区。这种方法需要更好的重新编码方法,对于K> 4是实用的。这绝对不是一种最好的或者优选的方法。然而,它在概念上简单并且容易被认为是从根本上无偏见的。它对于小K来说也非常快。脚本在Sage笔记本上运行良好,并且不会调用Sage函数。它不是随机抽样的脚本。随机抽样本身不是问题。

的方法:

1)治疗整数分区仿佛其被加数是根据在第一词汇的分区,例如最大被加数的大小用零串接在一起和填充[17,1,1,1] - > 17010101 & [5,5,5,5] - > 05050505

2.)将结果整数视为从最大整数中减去(即int代表第一个词法分区)。例如3.将每个得到的减小的整数除以一个公分母,例如, 99分之11959596= 120804

所以,如果我们要选择我们会随机分区:

1)选择0和120804(而非5050505和17010101)

之间的数字之间的数字

2.)从17010101

3. 99和相乘。减去数)根据我们如何处理每个整数为0的

被填充拆分得到的整数Pro's和Con's:正如问题的主体所述,这种特定的记录方法不足以大大提高随机选择代表P成员的整数的机会。例如K和实质上更大的总数,例如, N> 100,实现这个概念的函数可以非常快,因为该方法避免了及时递归(蛇吞食它的尾巴),这减慢了其他随机分区函数或使得其他函数对于处理大N不切实际。

在小K ,当考虑过程的其余部分有多快时,抽取P成员的概率可能是合理的。结合快速随机抽取,解码和评估,该功能可以找到用其他算法无法解决的N & K(例如N = 20000,K = 4)的组合的均匀随机分区。 一个更好的方法来重新编码整数是非常需要的,使这是一个通常强大的方法。

import random 
import sys 

首先,一些通常是有用的和简单的功能

def first_partition(N,K): 
    part = [N-K+1] 
    ones = [1]*(K-1) 
    part.extend(ones) 
return part 

def last_partition(N,K): 
    most_even = [int(floor(float(N)/float(K)))]*K 
    _remainder = int(N%K) 
    j = 0 

    while _remainder > 0: 
     most_even[j] += 1 
     _remainder -= 1 
     j += 1 
    return most_even 

def first_part_nmax(N,K,Nmax): 
    part = [Nmax] 
    N -= Nmax 
    K -= 1 
    while N > 0: 
     Nmax = min(Nmax,N-K+1) 
     part.append(Nmax) 
     N -= Nmax 
     K -= 1 

    return part 


#print first_partition(20,4) 
#print last_partition(20,4) 
#print first_part_nmax(20,4,12) 
#sys.exit() 

def portion(alist, indices): 
    return [alist[i:j] for i, j in zip([0]+indices, indices+[None])] 

def next_restricted_part(part,N,K): # *find next partition matching N&K w/out recursion 

    if part == last_partition(N,K):return first_partition(N,K) 
    for i in enumerate(reversed(part)): 
     if i[1] - part[-1] > 1: 
      if i[0] == (K-1): 
       return first_part_nmax(N,K,(i[1]-1)) 

      else: 
       parts = portion(part,[K-i[0]-1]) # split p 
       h1 = parts[0] 
       h2 = parts[1] 
       next = first_part_nmax(sum(h2),len(h2),(h2[0]-1)) 
       return h1+next 

""" *I don't know a math software that has this function and Nijenhuis and Wilf (1978) 
don't give it (i.e. NEXPAR is not restricted by K). Apparently, folks often get the 
next restricted part using recursion, which is unnecessary """ 

def int_to_list(i): # convert an int to a list w/out padding with 0' 
    return [int(x) for x in str(i)] 

def int_to_list_fill(i,fill):# convert an int to a list and pad with 0's 
    return [x for x in str(i).zfill(fill)] 

def list_to_int(l):# convert a list to an integer 
    return "".join(str(x) for x in l) 

def part_to_int(part,fill):# convert an int to a partition of K parts 
# and pad with the respective number of 0's 

    p_list = [] 
    for p in part: 
     if len(int_to_list(p)) != fill: 
      l = int_to_list_fill(p,fill) 
      p = list_to_int(l) 
     p_list.append(p) 
    _int = list_to_int(p_list) 
    return _int  

def int_to_part(num,fill,K): # convert an int to a partition of K parts 
    # and pad with the respective number of 0's 
    # This function isn't called by the script, but I thought I'd include 
    # it anyway because it would be used to recover the respective partition 

    _list = int_to_list(num) 
    if len(_list) != fill*K: 
     ct = fill*K - len(_list) 
     while ct > 0:  
      _list.insert(0,0) 
      ct -= 1  
    new_list1 = [] 
    new_list2 = [] 

    for i in _list: 
     new_list1.append(i) 
     if len(new_list1) == fill: 
      new_list2.append(new_list1) 
      new_list1 = [] 

    part = [] 
    for i in new_list2: 
     j = int(list_to_int(i)) 
     part.append(j) 

    return part 

最后,我们得到的全氮和零件数量K.下面将打印时,满足N & K的分区词汇顺序,带有相关的重新编码整数

N = 20 
K = 4 

print '#, partition, coded, _diff, smaller_diff' 

first_part = first_partition(N,K) # first lexical partition for N&K 
fill = len(int_to_list(max(first_part))) 
# pad with zeros to 1.) ensure a strictly decreasing relationship w/in P, 
# 2.) keep track of (encode/decode) partition summand values 

first_num = part_to_int(first_part,fill) 
last_part = last_partition(N,K) 
last_num = part_to_int(last_part,fill) 

print '1',first_part,first_num,'',0,'  ',0 

part = list(first_part) 
ct = 1 
while ct < 10: 
    part = next_restricted_part(part,N,K) 
    _num = part_to_int(part,fill) 

    _diff = int(first_num) - int(_num) 
    smaller_diff = (_diff/99) 
    ct+=1 

    print ct, part, _num,'',_diff,' ',smaller_diff 

OUTPUT:

CT,分区,编码,_diff,smaller_diff

1 [17,1,1,1] 17010101 0 0

2 [16,2,1,1] 16020101 990000 10000

3 [15,3,1,1] 15030101 1980000 20000

4 [15,2,2,1] 15020201 1989900 20100

5 [14,4,1,1] 14040101 2970000 30000

6 [14,3,2,1] 14030201 2979900 30100

7 [14,2,2,2] 14020202 2989899 30201

8 [13,5,1,1] 13050101 3960000 40000

9 [13,4,2,1] 13040201 3969900 40100

10 [13,3,3,1] 13030301 3979800 40200

I总之,最后一列的整数可能会小很多。

为什么基于这一想法的随机取样的策略是从根本上无偏:

N个具有K个每个份整数分区对应于一个且仅一个重新编码的整数。也就是说,我们不会随机挑选一个数字,对其进行解码,然后尝试重新排列元素以形成适当的N分区。因此,每个整数(无论是否对应于N & K的分区)都具有同样的机会被吸引。目标是固有地减少与K个部分不相对应的整数数量,从而使随机采样过程更快。

+1

我仍然有点不清楚为什么你觉得正确的方法是改善你的连接编码,而不是提出一个更好的初始编码。它看起来像你以前的[解决方案](http://stackoverflow.com/questions/10287021/an-algorithm-for-randomly-generating-integer-partitions-of-a-particular-length/12742508#12742508)以及我在这里提出的解决方案都计算N与K部分的分区之间的双射,[1..M]其中M是这样的分区的总数。你不能得到更紧凑的编码。 – mhum 2013-02-18 21:22:33

+0

我不会说我觉得这是正确的做法。这只是一种我用来生成N部分的随机分区的方法,其中K部分没有递归,这是简单的,从根本上无偏见的,而且没有偏见。当然,当K> 4时,这是疯狂的不切实际的。再一次,这个想法是递归是一些随机分割问题的敌人。 – klocey 2013-02-18 21:30:59

+0

您以前的解决方案和我在这里概述的解决方案基本相同,除了我使用稍微不同的递归。如果您按照我所做的操作并且只是迭代地预先计算D,那么您在执行以前的解决方案(速度,堆栈碎片)时遇到的问题可能会得到缓解。 – mhum 2013-02-18 21:41:20