logo标签列表

使用给定的（数字）分布生成随机数

pythonmodulerandom
236

我有一个包含不同值概率的文件，例如：

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

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

我需要这个是因为我想生成一组生日（它们不遵循标准random模块中的任何分布）。

2
除了random.choice()之外，您可以使用适当数量的出现次数构建主列表并选择一个。当然，这是一个重复的问题。 - S.Lott
2
@S.Lott对于分布差异较大的情况，这是否会占用大量内存？ - Lucas Moeskops
2
@S.Lott：对于少量出现次数，您的选择方法可能很好，但在不必要时我宁愿避免创建巨大的列表。 - pafcu
8
@S.Lott: 好的，大约有10000365=3650000=360万个元素。我不确定Python中的内存使用情况，但至少是3.6M4B=14.4MB。这不是很多，但在存在一个同样简单且不需要额外内存的方法时，你也不应该忽略它。 - pafcu
1
如果概率不是有理数，@S.Lott就无法工作。如果我有两个项目，一个概率为sqrt(2)/2，另一个概率为1-sqrt(2)/2，如果我想从这个分布中采样，相对频率需要达到10位小数精度，那么我需要有一个包含大约10^10个重复项的主表。有更高效的方法来做到这一点，而且它们非常简单，不需要主表。 - Rafael S. Calsaverini
显示剩余9条评论
13个回答
207

scipy.stats.rv_discrete可能是你想要的。你可以通过values参数提供概率值。然后，你可以使用分布对象的rvs()方法生成随机数。

正如Eugene Pakhomov在评论中指出的那样，你还可以通过p关键字参数传递给numpy.random.choice()，例如。

numpy.random.choice(numpy.arange(1, 7), p=[0.1, 0.05, 0.05, 0.2, 0.4, 0.2])

如果您使用的是Python 3.6或更高版本，您可以使用标准库中的random.choices()函数 - 请参考Mark Dickinson的回答。
18
在我的电脑上，numpy.random.choice() 函数快了将近20倍。 - Eugene Pakhomov
@EugenePakhomov，我不太理解你的评论。所以一个完全不同的函数比我建议的那个更快。我的建议仍然是使用能够实现你想要的功能的函数，而不是一个做其他事情但速度更快的函数。 - Sven Marnach
11
它在原始问题方面执行完全相同的操作。例如：numpy.random.choice(numpy.arange(1, 7), p=[0.1, 0.05, 0.05, 0.2, 0.4, 0.2]) - Eugene Pakhomov
2
令人惊讶的是，rv_discrete.rvs() 的运行时间和内存占用为 O(len(p) * size)！而 choice() 似乎在最优的 O(len(p) + log(len(p)) * size) 时间内运行。 - alyaxey
3
如果您正在使用Python 3.6或更高版本，则有另一个答案不需要任何附加包。 - Mark Ransom
194

自Python 3.6以来，Python标准库中有一个解决方案，即random.choices。

示例用法：让我们设置与OP问题中相匹配的种群和权重：

>>> from random import choices
>>> population = [1, 2, 3, 4, 5, 6]
>>> weights = [0.1, 0.05, 0.05, 0.2, 0.4, 0.2]

现在，choices(population, weights)生成一个包含单个样本的长度为1的列表：

>>> choices(population, weights)
[4]

可选的关键字参数k允许我们一次请求多个样本。这很有价值，因为每次调用random.choices之前都需要进行一些准备工作，以生成任何样本; 通过一次性生成多个样本，我们只需做一次准备工作。在这里，我们生成一百万个样本，并使用collections.Counter检查我们得到的分布大致符合我们给出的权重。

>>> million_samples = choices(population, weights, k=10**6)
>>> from collections import Counter
>>> Counter(million_samples)
Counter({5: 399616, 6: 200387, 4: 200117, 1: 99636, 3: 50219, 2: 50025})
这个有 Python 2.7 版本吗？ - abbas786
1
@abbas786：不是内置的，但是这个问题的其他答案都应该在Python 2.7上工作。如果愿意的话，您还可以查找Python 3的random.choices源代码并复制它。 - Mark Dickinson
1
对我来说，random.choices 与 k=1 结合使用会返回一个长度为一的列表，即 choices(population, weights) 应该返回 [4]。 - christianbrodbeck
@christianbrodbeck：谢谢，已修复。我几乎总是通过复制和粘贴来生成这些片段，所以显然这里出了问题。 - Mark Dickinson
谢谢！我一直在想这是否是版本问题，但这解释了一切。 - christianbrodbeck
有没有办法确保样本中不会重复包含同一个成员？ - theonlygusti
33

使用累积分布函数生成列表的一个优点是可以使用二分搜索。虽然需要 O(n) 的预处理时间和空间，但可以在 O(k log n) 的时间内获取 k 个数字。由于普通的 Python 列表效率低下，因此可以使用 array 模块。

如果您坚持要使用恒定空间，可以采取以下措施； 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
在你的实现中，(item, prob) 对的顺序在列表中很重要，对吗？ - stackoverflowuser2010
1
@stackoverflowuser2010：这不应该有影响（除了浮点数误差）。 - sdcvvc
很好。我发现这比scipy.stats.rv_discrete快30%。 - Adrienne
1
这个函数很多时候会因为最后一行而抛出 KeyError 的异常。 - imrek
@DrunkenMaster：我不明白。你知道 l[-1] 返回列表的最后一个元素吗？ - sdcvvc
显示剩余4条评论
18

好的，我知道你正在寻找最简化的解决方案，但也许那些自己制定的解决方案不够简洁，以符合你的喜好。 :-)

pdf = [(1, 0.1), (2, 0.05), (3, 0.05), (4, 0.2), (5, 0.4), (6, 0.2)]
cdf = [(i, sum(p for j,p in pdf if j < i)) for i,_ in pdf]
R = max(i for r in [random.random()] for i,c in cdf if c <= r)

我通过仔细观察这个表达式的输出，伪确认这个方法是可行的：

sorted(max(i for r in [random.random()] for i,c in cdf if c <= r)
       for _ in range(1000))
1
这看起来很不错。为了让事情更清晰，这里是上述代码连续执行3次的结果：['概率为0.1的1的数量为：113', '概率为0.05的2的数量为：55', '概率为0.05的3的数量为：50', '概率为0.2的4的数量为：201', '概率为0.4的5的数量为：388', '概率为0.2的6的数量为：193']..............['概率为0.1的1的数量为：77', '概率为0.05的2的数量为：60', '概率为0.05的3的数量为：51', '概率为0.2的4的数量为：193', '概率为0.4的5的数量为：438', '概率为0.2的6的数量为：181'] ............. - Vaibhav
1 的概率为 0.1 的数量为 84，2 的概率为 0.05 的数量为 52，3 的概率为 0.05 的数量为 53，4 的概率为 0.2 的数量为 210，5 的概率为 0.4 的数量为 405，6 的概率为 0.2 的数量为 196。 - Vaibhav
一个问题，如果“i”是一个对象，我该如何返回max(i)？ - Vaibhav
@Vaibhav i 不是一个对象。 - Marcelo Cantos
18
也许现在有点晚了，但您可以使用numpy.random.choice()函数，并传递p参数：
val = numpy.random.choice(numpy.arange(1, 7), p=[0.1, 0.05, 0.05, 0.2, 0.4, 0.2])
1
OP 不想使用 random.choice() - 请查看评论。 - pobrelkey
6
numpy.random.choice()与random.choice()截然不同，支持概率分布。 - Eugene Pakhomov
我不能使用函数来定义p吗？为什么要用数字来定义它呢？ - rjurney
如果您想从特定分布中进行抽样，应该使用像scipy.stats或statsmodels这样的统计软件包，然后从您想要抽样的具体概率分布中获取样本。此问题涉及用户定义的离散分布情况。 - Heberto Mayorquin
10

我为从自定义连续分布中绘制随机样本编写了一种解决方案。

我需要这个解决方案来处理类似于您的用例（即使用给定概率分布生成随机日期）。

您只需要函数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()

连续定制分布和离散样本分布

这个解决方案的性能肯定可以提高，但我更倾向于可读性。

为什么连续分布的密度会小于1？ - erwan gaymard
2

根据它们的权重列出物品清单：

items = [1, 2, 3, 4, 5, 6]
probabilities= [0.1, 0.05, 0.05, 0.2, 0.4, 0.2]
# if the list of probs is normalized (sum(probs) == 1), omit this part
prob = sum(probabilities) # find sum of probs, to normalize them
c = (1.0)/prob # a multiplier to make a list of normalized probs
probabilities = map(lambda x: c*x, probabilities)
print probabilities

ml = max(probabilities, key=lambda x: len(str(x)) - str(x).find('.'))
ml = len(str(ml)) - str(ml).find('.') -1
amounts = [ int(x*(10**ml)) for x in probabilities]
itemsList = list()
for i in range(0, len(items)): # iterate through original items
  itemsList += items[i:i+1]*amounts[i]

# choose from itemsList randomly
print itemsList

优化策略可能是通过最大公因数将金额标准化，以使目标列表更小。

此外，这篇文章可能会很有趣。

如果项目列表很大，这可能会使用大量额外的内存。 - pafcu
@pafcu 同意。只是一个解决方案，第二个浮现在我脑海中的（第一个是搜索类似“Python权重概率”的东西：））。 - khachik
1
from __future__ import division
import random
from collections import Counter


def num_gen(num_probs):
    # calculate minimum probability to normalize
    min_prob = min(prob for num, prob in num_probs)
    lst = []
    for num, prob in num_probs:
        # keep appending num to lst, proportional to its probability in the distribution
        for _ in range(int(prob/min_prob)):
            lst.append(num)
    # all elems in lst occur proportional to their distribution probablities
    while True:
        # pick a random index from lst
        ind = random.randint(0, len(lst)-1)
        yield lst[ind]

验证：
gen = num_gen([(1, 0.1),
               (2, 0.05),
               (3, 0.05),
               (4, 0.2),
               (5, 0.4),
               (6, 0.2)])
lst = []
times = 10000
for _ in range(times):
    lst.append(next(gen))
# Verify the created distribution:
for item, count in Counter(lst).iteritems():
    print '%d has %f probability' % (item, count/times)

1 has 0.099737 probability
2 has 0.050022 probability
3 has 0.049996 probability 
4 has 0.200154 probability
5 has 0.399791 probability
6 has 0.200300 probability
1

基于其他解决方案，您可以生成累积分布（无论是整数还是浮点数），然后可以使用二分查找使其更快。

这是一个简单的例子（我在这里使用了整数）。

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

现在我们使用random.randint随机选择一个数字，范围在20+60+10+10以内，然后使用bisect快速获取实际值

1

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

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
distribution列表需要按概率排序吗？ - YQ.Wang
如果按照概率从大到小排序，它不一定需要，但是它会执行得最快。 - Lucas Moeskops
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接