Python：从几何分布生成数据

Question

这是生成随机数的最佳或最有效的方法吗？这些随机数来自一个几何分布，并且参数数组中可能包含0。

allids["c"]=[2,0,1,1,3,0,0,2,0]
[ 0 if x == 0 else numpy.random.geometric(1./x) for x in allids["c"]]

我有点担心优化的问题。

编辑：

稍微解释一下背景：我有一串字符（比如 ATCGGGA），我想要扩展或缩短某个字符的连续出现次数（比如，如果原始序列中有两个 'A' 连续出现，我想模拟一个序列，使得这个序列中 'A' 的期望值也是2，但根据几何分布来变化）。而那些只出现一次的字符我不想让它们的长度变化。

所以如果

seq = 'AATCGGGAA'
allids["c"]=[2,0,1,1,3,0,0,2,0]
rep=[ 0 if x == 0 else numpy.random.geometric(1./x) for x in allids["c"]]

"".join([s*r for r, s in zip(rep, seq)])

当 rep 是 [1, 0, 1, 1, 3, 0, 0, 1, 0] 时，将会输出

"ATCGGGA"

Answer 1

这个怎么样：

counts = array([2, 0, 1, 1, 3, 0, 0, 2, 0], dtype=float)
counts_ma = numpy.ma.array(counts, mask=(counts == 0))
counts[logical_not(counts.mask)] = \
    array([numpy.random.geometric(v) for v in 1.0 / counts[logical_not(counts.mask)]])

你可以提前计算好同种聚合物的分布情况，这样就可以减少对geometric的调用次数。因为从随机数生成器一次性获取大量值，比每次单独调用要高效得多。

Answer 2

你可以使用一个叫做 掩码数组 的东西来避免除以零的问题。

import numpy as np
a = np.ma.masked_equal([2, 0, 1, 1, 3, 0, 0, 2, 0], 0)
rep = np.random.geometric(1. / a)
rep[a.mask] = 0

这个方法会为每个a 的元素生成一个随机样本，然后再把其中一些删除。如果你担心这样会浪费随机数，你可以只生成刚好够用的数量，像这样：

import numpy as np
a = np.ma.masked_equal([2, 0, 1, 1, 3, 0, 0, 2, 0], 0)
rep = np.zeros(a.shape, dtype=int)
rep[~a.mask] = np.random.geometric(1. / a[~a.mask])