2016-03-01 56 views
4

我想在numpy中运行一个相对简单的随机绘制,但我找不到表达它的好方法。 我认为最好的方法是将其描述为从没有更换的瓮中绘制。我有一个有k种颜色的瓮和每种颜色的n_k球。我想画出米球,并知道我有多少种颜色的球。来自瓮的Numpy绘图

我的当前的尝试它

np.bincount(np.random.permutation(np.repeat(np.arange(k), n_k))[:m], minlength=k)

这里,n_k是长度为k的阵列与所述球的计数。

看来,这相当于 np.bincount(np.random.choice(k, m, n_k/n_k.sum(), minlength=k)

这是一个好一点,但仍然不是很大。

回答

8

你想要什么multivariate hypergeometric distribution的实现。我不知道一个在numpy或scipy中,但它可能已经存在了某处。

您可以使用重复调用numpy.random.hypergeometric来实现它。这是否会比实施效率更高取决于有多少种颜色以及每种颜色有多少个球。

例如,下面的打印含有三种颜色瓮抽选结果的脚本(红,绿,蓝):

from __future__ import print_function 

import numpy as np 


nred = 12 
ngreen = 4 
nblue = 18 

m = 15 

red = np.random.hypergeometric(nred, ngreen + nblue, m) 
green = np.random.hypergeometric(ngreen, nblue, m - red) 
blue = m - (red + green) 

print("red: %2i" % red) 
print("green: %2i" % green) 
print("blue: %2i" % blue) 

输出示例:

red: 6 
green: 1 
blue: 8 

下面的函数概括地说,选择m球给出了一个数组colors保存着每种颜色的编号:

def sample(m, colors): 
    """ 
    Parameters 
    ---------- 
    m : number balls to draw from the urn 
    colors : one-dimensional array of number balls of each color in the urn 

    Returns 
    ------- 
    One-dimensional array with the same length as `colors` containing the 
    number of balls of each color in a random sample. 
    """ 

    remaining = np.cumsum(colors[::-1])[::-1] 
    result = np.zeros(len(colors), dtype=np.int) 
    for i in range(len(colors)-1): 
     if m < 1: 
      break 
     result[i] = np.random.hypergeometric(colors[i], remaining[i+1], m) 
     m -= result[i] 
    result[-1] = m 
    return result 

例如,

>>> sample(10, [2, 4, 8, 16]) 
array([2, 3, 1, 4]) 
+0

很酷。非常有趣的解决方案,谢谢如果只有多项式超几何;) –

+0

对,“样本”基本上是多变量超几何分布的实现(https://en.wikipedia.org/wiki/Hypergeometric_distribution#Multivariate_hypergeometric_distribution)。 –

+1

请注意,pymc提供了一个函数('rmultivariate_hypergeometric')用于从多元超几何绘制随机变量。链接在这里:http://pymc-devs.github.io/pymc/_modules/pymc/distributions.html – jvd10

0

下面应该工作:

def make_sampling_arr(n_k): 
    out = [ x for s in [ [i] * n_k[i] for i in range(len(n_k)) ] for x in s ] 
    return out 

np.random.choice(make_sampling_arr(n_k), m) 
+0

我要计数,然后你需要在最后调用bincount。 ''make_sampling_arr''和'np.repeat''一样,就像我上面所用的那样,对吧? –