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