用给定的(数字)分布生成随机数字
我有一个文件,有不同的值的概率,例如:
1 0.1 2 0.05 3 0.05 4 0.2 5 0.4 6 0.2
我想用这个分布生成随机数。 处理这个的现有模块是否存在? 编写自己的代码非常简单(构build累积密度函数,生成一个随机值[0,1]并select相应的值),但似乎这应该是一个常见问题,可能有人创build了一个函数/模块它。
我需要这个,因为我想生成一个生日的列表(不遵循标准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])
使用CDF生成列表的优点是可以使用二分search。 当需要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
(好吧,我知道你们正在寻求收缩包装,但是也许这些自制的解决scheme并不足以满足你的喜好。:-)
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)
我通过观察这个expression式的输出来伪证实了这一点:
sorted(max(i for r in [random.random()] for i,c in cdf if c <= r) for _ in range(1000))
也许这是晚了。 但是你可以使用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])
自Python 3.6以来,在Python的标准库中有一个解决scheme,即random.choices
。
示例用法:让我们build立一个匹配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)
生成一个样本:
>>> 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})
你可能想看看NumPy 随机抽样分布
根据weights
制定项目清单:
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
优化可以是通过最大公约数来标准化量,以使目标列表更小。
另外, 这可能是有趣的。
另一个答案,可能更快:)
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
基于其他解决scheme,你可以产生累积分布(如你喜欢的整数或浮点数),那么你可以使用平分来快速
这是一个简单的例子(我在这里使用整数)
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
函数将它从get_cdf
转换为20,20 + get_cdf
+ 60 + 10,20 + 60 + 10 + 10
现在我们用random.randint
select一个20 + 60 + 10 + 10的随机数,然后我们用对分来快速得到实际值
这些答案都不是特别清楚或简单。
这是一个明确,简单的方法,保证工作。
accumulate_normalize_probabilities需要一个将符号映射到概率或频率的字典p
。 它输出可用于select的元组列表的可用列表。
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)。 这些时间间隔可用于从列表中简单地逐步select(并从而对所提供的分布进行采样),直到间隔0.0 – > 1.0(先前准备)的随机数小于或等于当前符号的间隔终点。
规范化使我们摆脱了确保一切都达到一定价值的需要。 归一化之后,概率的“向量”总和为1.0。
下面的代码供select和生成任意长度的样本:
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
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]
validation:
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
这是一个更有效的方法 :
只需要调用你的'权重'数组(假设索引作为相应的项目)和没有。 需要的样品。 这个function可以很容易地修改来处理有序对。
返回索引(或项目)采样/采摘(与replace)使用各自的概率:
def resample(weights, n): beta = 0 # Caveat: Assign max weight to max*2 for best results max_w = max(weights)*2 # Pick an item uniformly at random, to start with current_item = random.randint(0,n-1) result = [] for i in range(n): beta += random.uniform(0,max_w) while weights[current_item] < beta: beta -= weights[current_item] current_item = (current_item + 1) % n # cyclic else: result.append(current_item) return result
关于while循环中使用的概念的简短说明。 我们从累积贝塔(这是一个随机均匀构造的累积值)中减去当前物品的重量,并增加当前索引以find与贝塔值相匹配的物品。