我有一个文件,不同的值的一些概率,例如:

1 0.1
2 0.05
3 0.05
4 0.2
5 0.4
6 0.2

我想用这个分布生成随机数。是否存在处理此问题的现有模块?自己编写代码是相当简单的(构建累积密度函数,生成一个随机值[0,1]并选择相应的值),但这似乎应该是一个常见的问题,可能有人已经为它创建了一个函数/模块。

我需要这个,因为我想生成一个生日列表(它不遵循标准随机模块中的任何分布)。


当前回答

使用CDF生成列表的一个优点是可以使用二分搜索。当你需要O(n)个时间和空间进行预处理时,你可以得到O(k log n)个数字。由于普通的Python列表效率很低,你可以使用数组模块。

如果你坚持空间不变,你可以做到以下几点;O(n)时间,O(1)空间。

def random_distr(l):
    r = random.uniform(0, 1)
    s = 0
    for item, prob in l:
        s += prob
        if s >= r:
            return item
    return item  # Might occur because of floating point inaccuracies

其他回答

基于其他解决方案,您可以生成累积分布(作为整数或浮点数),然后您可以使用平分使其更快

这是一个简单的例子(我在这里使用整数)

l=[(20, 'foo'), (60, 'banana'), (10, 'monkey'), (10, 'monkey2')]
def get_cdf(l):
    ret=[]
    c=0
    for i in l: c+=i[0]; ret.append((c, i[1]))
    return ret

def get_random_item(cdf):
    return cdf[bisect.bisect_left(cdf, (random.randint(0, cdf[-1][0]),))][1]

cdf=get_cdf(l)
for i in range(100): print get_random_item(cdf),

get_cdf函数会将20、60、10、10转换为20、20+60、20+60+10、20+60+10+10

现在我们随机选择一个20+60+10+10的随机数。然后我们用二分法快速得到实际值

你可能想看看NumPy随机抽样分布

另一个答案,可能更快:)

distribution = [(1, 0.2), (2, 0.3), (3, 0.5)]  
# init distribution  
dlist = []  
sumchance = 0  
for value, chance in distribution:  
    sumchance += chance  
    dlist.append((value, sumchance))  
assert sumchance == 1.0 # not good assert because of float equality  

# get random value  
r = random.random()  
# for small distributions use lineair search  
if len(distribution) < 64: # don't know exact speed limit  
    for value, sumchance in dlist:  
        if r < sumchance:  
            return value  
else:  
    # else (not implemented) binary search algorithm  

这些答案都不是特别明确或简单的。

这里有一个明确、简单、保证有效的方法。

accumulate_normalize_probability接受一个字典p,将符号映射到概率或频率。它输出可用的元组列表,从中进行选择。

def accumulate_normalize_values(p):
        pi = p.items() if isinstance(p,dict) else p
        accum_pi = []
        accum = 0
        for i in pi:
                accum_pi.append((i[0],i[1]+accum))
                accum += i[1]
        if accum == 0:
                raise Exception( "You are about to explode the universe. Continue ? Y/N " )
        normed_a = []
        for a in accum_pi:
                normed_a.append((a[0],a[1]*1.0/accum))
        return normed_a

收益率:

>>> accumulate_normalize_values( { 'a': 100, 'b' : 300, 'c' : 400, 'd' : 200  } )
[('a', 0.1), ('c', 0.5), ('b', 0.8), ('d', 1.0)]

为什么它有效

累积步骤将每个符号转换为它自身与前一个符号的概率或频率之间的间隔(或第一个符号的情况为0)。这些间隔可以通过简单地遍历列表,直到间隔0.0 -> 1.0(前面准备的)中的随机数小于或等于当前符号的间隔终点来进行选择(从而对所提供的分布进行抽样)。

规范化使我们不再需要确保所有内容的总和为某个值。归一化后,概率的“向量”总和为1.0。

从分布中选择和生成任意长样本的其余代码如下:

def select(symbol_intervals,random):
        print symbol_intervals,random
        i = 0
        while random > symbol_intervals[i][1]:
                i += 1
                if i >= len(symbol_intervals):
                        raise Exception( "What did you DO to that poor list?" )
        return symbol_intervals[i][0]


def gen_random(alphabet,length,probabilities=None):
        from random import random
        from itertools import repeat
        if probabilities is None:
                probabilities = dict(zip(alphabet,repeat(1.0)))
        elif len(probabilities) > 0 and isinstance(probabilities[0],(int,long,float)):
                probabilities = dict(zip(alphabet,probabilities)) #ordered
        usable_probabilities = accumulate_normalize_values(probabilities)
        gen = []
        while len(gen) < length:
                gen.append(select(usable_probabilities,random()))
        return gen

用法:

>>> gen_random (['a','b','c','d'],10,[100,300,400,200])
['d', 'b', 'b', 'a', 'c', 'c', 'b', 'c', 'c', 'c']   #<--- some of the time

我写了一个从自定义连续分布中抽取随机样本的解决方案。

我需要这个类似于你的用例(即生成随机日期与给定的概率分布)。

你只需要函数random_custDist和行samples=random_custDist(x0,x1,custDist=custDist,size=1000)。其余的都是装饰^^。

import numpy as np

#funtion
def random_custDist(x0,x1,custDist,size=None, nControl=10**6):
    #genearte a list of size random samples, obeying the distribution custDist
    #suggests random samples between x0 and x1 and accepts the suggestion with probability custDist(x)
    #custDist noes not need to be normalized. Add this condition to increase performance. 
    #Best performance for max_{x in [x0,x1]} custDist(x) = 1
    samples=[]
    nLoop=0
    while len(samples)<size and nLoop<nControl:
        x=np.random.uniform(low=x0,high=x1)
        prop=custDist(x)
        assert prop>=0 and prop<=1
        if np.random.uniform(low=0,high=1) <=prop:
            samples += [x]
        nLoop+=1
    return samples

#call
x0=2007
x1=2019
def custDist(x):
    if x<2010:
        return .3
    else:
        return (np.exp(x-2008)-1)/(np.exp(2019-2007)-1)
samples=random_custDist(x0,x1,custDist=custDist,size=1000)
print(samples)

#plot
import matplotlib.pyplot as plt
#hist
bins=np.linspace(x0,x1,int(x1-x0+1))
hist=np.histogram(samples, bins )[0]
hist=hist/np.sum(hist)
plt.bar( (bins[:-1]+bins[1:])/2, hist, width=.96, label='sample distribution')
#dist
grid=np.linspace(x0,x1,100)
discCustDist=np.array([custDist(x) for x in grid]) #distrete version
discCustDist*=1/(grid[1]-grid[0])/np.sum(discCustDist)
plt.plot(grid,discCustDist,label='custom distribustion (custDist)', color='C1', linewidth=4)
#decoration
plt.legend(loc=3,bbox_to_anchor=(1,0))
plt.show()

这个解决方案的性能肯定是可以改进的,但我更喜欢可读性。