有效使用 numpy.random.choice 与重复数字和替代方案

Question

我需要生成一个包含重复元素的大数组，我的代码是：

np.repeat(xrange(x,y), data)

但是，data 是一个类型为 float64 的 numpy 数组（但它代表整数，那里没有 2.1），我得到了错误

TypeError: Cannot cast array data from dtype('float64') to dtype('int64') according to the rule 'safe'

示例：

In [35]: x
Out[35]: 26

In [36]: y
Out[36]: 50

In [37]: data
Out[37]: 
array([ 3269.,   106.,  5533.,   317.,  1512.,   208.,   502.,   919.,
     406.,   421.,  1690.,  2236.,   705.,   505.,   230.,   213.,
     307.,  1628.,  4389.,  1491.,   355.,   103.,   854.,   424.])
In [38]: np.repeat(xrange(x,y), data)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call    last)
<ipython-input-38-105860821359> in <module>()
----> 1 np.repeat(xrange(x,y), data)

/home/pcadmin/anaconda2/lib/python2.7/site-packages/numpy    /core/fromnumeric.pyc in repeat(a, repeats, axis)
394         repeat = a.repeat
395     except AttributeError:
--> 396         return _wrapit(a, 'repeat', repeats, axis)
397     return repeat(repeats, axis)
398 

/home/pcadmin/anaconda2/lib/python2.7/site-packages/numpy  /core/fromnumeric.pyc in _wrapit(obj, method, *args, **kwds)
 46     except AttributeError:
 47         wrap = None
---> 48     result = getattr(asarray(obj), method)(*args, **kwds)
 49     if wrap:
 50         if not isinstance(result, mu.ndarray):

TypeError: Cannot cast array data from dtype('float64') to dtype('int64') according to the rule 'safe'

我通过将代码更改为来解决它

np.repeat(xrange(x,y), data.astype('int64'))

但是，这现在是我的代码中最昂贵的行之一！还有其他选择吗？

顺便说一下，我在里面用这个

np.random.choice(np.repeat(xrange(x,y), data.astype('int64')), z)

为了得到一个样本而不用替换 x 和 y 之间的整数的大小 z，每个样本的数量在数据中给出。我想这也是最好的方法吧？

Answer 1

潜伏在问题中的是multivariate hypergeometric distribution。在Numpy drawing from urn 中，我实现了一个从这个分布中抽取样本的函数。我怀疑它与答案中描述的解决方案@DiogoSantos 非常相似。 Diogo 说使用这种方法很慢，但我发现以下方法比 Divakar 的 optmized_v1 更快。

这是一个函数，它使用链接答案中的sample(n, colors) 来实现与 Divakar 函数具有相同签名的函数。

def hypergeom_version(x, y, z, data):
    s = sample(z, data)
    result = np.repeat(np.arange(x, y), s)
    return result

（这会以排序的顺序返回值。如果您需要将值按随机顺序排列，请在返回语句之前添加np.random.shuffle(result)。它不会显着改变执行时间。）

比较：

In [153]: x = 100

In [154]: y = 100100

In [155]: z = 10000

In [156]: data = np.random.randint(1, 125, (y-x)).astype(float)

Divakar 的optimized_v1：

In [157]: %timeit optimized_v1(x, y, z, data)
1 loop, best of 3: 520 ms per loop

hypergeom_version:

In [158]: %timeit hypergeom_version(x, y, z, data)
1 loop, best of 3: 244 ms per loop

data中的值越大，相对性能越好：

In [164]: data = np.random.randint(100, 500, (y-x)).astype(float)

In [165]: %timeit optimized_v1(x, y, z, data)
1 loop, best of 3: 2.91 s per loop

In [166]: %timeit hypergeom_version(x, y, z, data)
1 loop, best of 3: 246 ms per loop

Answer 2

问题陈述

这个问题很有趣！只是为了让读者了解问题而不涉及次要的数据转换问题，我们有一系列值，比如说a = np.arange(5)，即

a = np.array([0,1,2,3,4])

现在，假设我们有另一个数组，其中列出了 a 中每个 5 数字的重复次数。所以，让他们成为：

reps = np.array([2,4,6,2,2])

接下来，我们正在执行这些重复：

In [32]: rep_nums = np.repeat(a,reps)

In [33]: rep_nums
Out[33]: array([0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 4, 4])

最后，我们希望使用np.random.choice() 从那些重复的数字中选择z 的元素数量，并且不进行替换。

假设z = 7 选择7 元素，那么使用np.random.choice()，我们将有：

In [34]: np.random.choice(rep_nums,7,replace=False)
Out[34]: array([2, 4, 0, 2, 4, 1, 2])

现在，这里的 without replacement 术语可能听起来令人困惑，因为我们已经在 rep_nums 中有重复的数字。但是，它本质上的意思是，np.random.choice() 的输出不能包含例如不止两个4's，因为rep_nums 有两个4's。

所以，问题是我们想要摆脱 np.repeat 部分，这可能是非常大的数组的瓶颈。

建议的方法

查看rep_nums 的输出，一种想法是生成z = 7 的唯一元素，其长度范围为rep_nums：

In [44]: np.random.choice(rep_nums.size,7,replace=False)
Out[44]: array([ 7,  2,  4, 10, 13,  8,  3])

这些数字表示该长度的索引。因此，我们只需要在rep_nums 中寻找每个7 数字将进入的bin（在5 bin 之外）。为此，我们可以使用np.searchsorted。因此，我们将有一个实现来处理通用的x、y，就像这样 -

# Get the intervals of those bins
intervals = data.astype(int).cumsum()

# Decide length of array if we had repeated with `np.repeat`
max_num = intervals[-1]

# Get unique numbers (indices in this case)
ids = np.random.choice(max_num,z,replace=False)

# Use searchsorted to get bin IDs and add in `x` offset
out = x+np.searchsorted(intervals,ids,'right')

运行时测试

功能：

def org_app(x,y,z,data):
    rep_nums = np.repeat(range(x,y), data.astype('int64'))
    out = np.random.choice(rep_nums, z,replace=False)
    return out
     
def optimized_v1(x,y,z,data):     
    intervals = data.astype(int).cumsum()
    max_num = intervals[-1]
    ids = np.random.choice(max_num,z,replace=False)
    out = x+np.searchsorted(intervals,ids,'right')
    return out

完整功能的时间安排 -

In [79]: # Setup inputs
    ...: x = 100
    ...: y = 10010
    ...: z = 1000
    ...: data = np.random.randint(100,5000,(y-x)).astype(float)
    ...: 

In [80]: %timeit org_app(x,y,z,data)
1 loop, best of 3: 7.17 s per loop

In [81]: %timeit optimized_v1(x,y,z,data)
1 loop, best of 3: 6.92 s per loop

看起来我们没有得到很好的加速。让我们深入挖掘一下，看看替换 np.repeat 节省了多少！

首先采用原始方法 -

In [82]: %timeit np.repeat(range(x,y), data.astype('int64'))
1 loop, best of 3: 227 ms per loop

让我们看看建议的方法在这方面取得了多大的改进。所以，让我们在提议的方法中为除np.random.choice() 之外的所有内容计时 -

In [83]: intervals = data.astype(int).cumsum()
    ...: max_num = intervals[-1]
    ...: ids = np.random.choice(max_num,z,replace=False)
    ...: out = x+np.searchsorted(intervals,ids,'right')
    ...: 

In [84]: %timeit data.astype(int).cumsum()
10000 loops, best of 3: 36.6 µs per loop

In [85]: %timeit intervals[-1]
10000000 loops, best of 3: 142 ns per loop

In [86]: %timeit x+np.searchsorted(intervals,ids,'right')
10000 loops, best of 3: 127 µs per loop

这比np.repeat的227ms好多了！！

因此，我们希望在非常大的阵列中，删除 np.repeat 的好处会真正发挥作用，否则 np.random.choice() 本身看起来就像是瓶颈。

Answer 3

为了完成，我还有一个替代实现。鉴于我们有data，我们可以对每个类使用超几何采样：

• 反向计算data.cumsum()
• 为每个班级抽奖np.hypergeometric(data[pos], cumsum[pos]-data[pos], remain)

但是，当我们有很多类且每个类中的单元很少时，这需要很长时间。