我正在为基因组学课程编写一些代码,但在某些方面遇到困难.
Hi I'm doing some code for a genomics class and I am having difficulty on a certain part.
我有一组相互排斥的事件 具有概率
I have a set of mutually exclusive events with probabilities
我想模拟以给定的概率对事件进行随机采样n次.
I want to simulate randomly sampling an event n times with the given probability.
输入:概率= {0.3,0.2,0.5}事件{e1,e2,e3} n = 100
input: probabilities = {0.3, 0.2, 0.5} events{e1,e2,e3} n=100
输出:e3应该有〜50个结果,e2应该有〜20个结果,而e1应该有〜30个结果. 请注意,这些可能不完全是50、20、30,因为 经验值不同于理论值...
output: there should be ~50 results for e3, ~20 for e2 and ~30 for e1. Note that these are probably not exactly 50, 20, 30 because empirical values are different from theoretical values...
推荐答案Python没有内置任何加权采样功能(NumPy/SciPy具备),但是对于像这样的非常简单的情况,这很简单:
Python doesn't have any weighted sampling functionality built in (NumPy/SciPy does), but for a really simple case like this, it's pretty easy:
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函数.如果清单真的很短,您可以用低效率的单线伪造它:
If you don't have Python 3.2+, you don't have the accumulate function; you can fake it with an inefficient one-liner if the list really is this short:
totals = [sum(probabilities[:i+1]) for i in range(len(probabilities))]…,或者您可以编写一个显式循环或丑陋的reduce调用,或从文档.
… or you can write an explicit loop, or an ugly reduce call, or copy the equivalent Python function from the docs.
此外,请注意,如果您可以确定数字加起来为1.0,则random.uniform(0, totals[-1])只是写random.random()的一种更为复杂的方式.
Also, note that random.uniform(0, totals[-1]) is just a more complicated way of writing random.random() if you can be sure that your numbers add up to 1.0.
一种快速的测试方法:
>>> samples = [sample() for _ in range(100000)] >>> samples.count(0) 29878 >>> samples.count(1) 19908 >>> samples.count(2) 50214分别接近100000的30%,20%和50%.
Those are pretty close to 30%, 20%, and 50% of 100000, respectively.