我有一个包含不同值概率的文件,例如:
1 0.1
2 0.05
3 0.05
4 0.2
5 0.4
6 0.2
我想使用这个分布生成随机数。是否存在处理此问题的现有模块?自己编写代码相当简单(构建累积密度函数,生成[0,1]的随机值并选择相应的值),但似乎这应该是一个常见问题,可能有人已经为此创建了函数/模块。
我需要这个是因为我想生成一组生日(它们不遵循标准random
模块中的任何分布)。
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])
random.choices()
函数 - 请参考Mark Dickinson的回答。numpy.random.choice()
函数快了将近20倍。 - Eugene Pakhomovnumpy.random.choice(numpy.arange(1, 7), p=[0.1, 0.05, 0.05, 0.2, 0.4, 0.2])
- Eugene Pakhomov自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})
random.choices
与 k=1
结合使用会返回一个长度为一的列表,即 choices(population, weights)
应该返回 [4]
。 - christianbrodbeck使用累积分布函数生成列表的一个优点是可以使用二分搜索。虽然需要 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
l[-1]
返回列表的最后一个元素吗? - sdcvvc好的,我知道你正在寻找最简化的解决方案,但也许那些自己制定的解决方案不够简洁,以符合你的喜好。 :-)
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))
i
不是一个对象。 - Marcelo Cantosnumpy.random.choice()
函数,并传递p
参数:val = numpy.random.choice(numpy.arange(1, 7), p=[0.1, 0.05, 0.05, 0.2, 0.4, 0.2])
random.choice()
- 请查看评论。 - pobrelkeynumpy.random.choice()
与random.choice()
截然不同,支持概率分布。 - Eugene Pakhomovscipy.stats
或statsmodels
这样的统计软件包,然后从您想要抽样的具体概率分布中获取样本。此问题涉及用户定义的离散分布情况。 - Heberto Mayorquin我为从自定义连续分布中绘制随机样本编写了一种解决方案。
我需要这个解决方案来处理类似于您的用例(即使用给定概率分布生成随机日期)。
您只需要函数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根据它们的权重
列出物品清单:
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
优化策略可能是通过最大公因数将金额标准化,以使目标列表更小。
此外,这篇文章可能会很有趣。
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
基于其他解决方案,您可以生成累积分布(无论是整数还是浮点数),然后可以使用二分查找使其更快。
这是一个简单的例子(我在这里使用了整数)。
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快速获取实际值
另一个答案,可能更快 :)
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
random.choice()
之外,您可以使用适当数量的出现次数构建主列表并选择一个。当然,这是一个重复的问题。 - S.Lott