嗨,我正在为基因组学课程编写一些代码,但在某些方面遇到了困难。
我有一组互斥事件 与概率
我想以给定的概率模拟随机采样n次事件。
输入:概率= {0.3,0.2,0.5}事件{e1,e2,e3} n = 100
输出:e3应该有〜50个结果,e2应该有〜20个结果,而e1应该有〜30个结果。请注意,这些可能不完全是50、20、30,因为经验值与理论值不同…
Python没有内置任何加权采样功能(NumPy / SciPy具有),但是对于这样一个非常简单的情况,这很简单:
import itertools import random probabilities = [0.3, 0.2, 0.5] totals = list(itertools.accumulate(probabilities)) def sample(): n = random.uniform(0, totals[-1]) for i, total in enumerate(totals): if n <= total: return i
如果您没有Python 3.2+,则没有此accumulate功能。如果清单确实很短,则可以使用低效率的单线伪造它:
accumulate
totals = [sum(probabilities[:i+1]) for i in range(len(probabilities))]
…,或者您可以编写一个显式循环或丑陋的reduce调用,或从docs复制等效的Python函数。
reduce
另外,请注意,如果可以确定数字加起来为1.0 ,random.uniform(0, totals[-1])则这只是一种更复杂的书写方式random.random()。
random.uniform(0, totals[-1])
random.random()
一种快速的测试方法:
>>> samples = [sample() for _ in range(100000)] >>> samples.count(0) 29878 >>> samples.count(1) 19908 >>> samples.count(2) 50214
这些分别分别接近100000的30%,20%和50%。