使用 Numpy 进行唯一随机数采样答案

【问题标题】：Unique random number sampling with Numpy使用 Numpy 进行唯一随机数采样
【发布时间】：2017-07-12 15:46:00
【问题描述】：

我需要创建一个 10,000 x 50 的数组，其中每一行包含 1 到 365 之间的一系列随机数，如下所示：

[[  4  11  14 ..., 355 360 364]
 [  2  13  15 ..., 356 361 361]
 [  4  12  18 ..., 356 361 365]
 ..., 
 [  6   9  17 ..., 356 362 364]
 [  1  10  19 ..., 352 357 360]
 [  1   9  17 ..., 356 358 364]]

我想出的唯一方法是通过迭代器：

sample_dates = np.array([np.sort(np.random.choice(365, 50, replace=False)) for _ in range(10000)])

这行得通，但速度很慢（运行约 0.33 秒），我将这样做数千次）。有没有更快的方法来做到这一点？

编辑：据我所知，这个解决方案最昂贵的部分是迭代和对 np.random.choice 的 10k 次单独调用，而不是排序

【问题讨论】：

你打算用这个做什么？
我建议您使用更快的方法，但我不想添加它作为答案，因为我不知道编程语言：在开头生成第一个随机数怎么样行，然后为第 n 个元素生成另一个随机数并将其添加到前一个随机数。你会得到一行排序的随机数。为了满足您的间隔要求，请使用 1 到 365/50 之间的随机数

标签： python performance numpy

【解决方案1】：

以下解决方案不使用排序：

l = np.array([True]*50 + [False]*315)
total = np.arange(1,366)
sample_dates = np.array([total[np.random.permutation(l)] for _ in range(10000)])

因此它似乎比其他建议的解决方案更快（在我的计算机上需要 0.44 秒，而“Nils Werner”的解决方案需要 0.77 秒。OP 的解决方案需要 0.81 秒）。

【讨论】：

【解决方案2】：

考虑到数组的形状，我认为对列进行迭代可能会带来一些改进。所以我的想法是生成 10k 个数字 - 替换。然后，在一个循环中，生成另外 10k 个数字并检查逐行重复。如果有的话，消除那些并生成那么多随机数。如果我没记错的话，这也称为命中和未命中算法。

这是工作代码：

arr = np.random.choice(365, 10000)
for i in range(49):
    arr2 = np.random.choice(365, 10000)
    comp = (arr2 == arr)
    while comp.any():
        duplicate = comp if i==0 else comp.any(axis=0)
        arr2[duplicate] = np.random.choice(365, duplicate.sum())
        comp = (arr2 == arr)
    arr = np.vstack([arr, arr2])
arr = arr.T
arr.sort(axis=1)

这需要 93.4 毫秒才能完成。您的尝试在我的计算机上花费了 590 毫秒，因此它提供了大约 6 倍的改进。

【讨论】：

【解决方案3】：

这是np.argpartition/np.argsort 的一个技巧。

这个想法是-

获取形状为(10000,365) 的随机浮点数组，沿每一行对其执行argsort。这将为我们提供唯一的索引，从而模拟与np.random.choice 一起使用的replace=False 标准。
切出每行的前50 列。
最后，对每一行进行排序来完成排序数据的工作。

现在，我们可以使用np.argpartition 进一步提高性能，以使用k=50 沿每一行进行分区。

因此，我们将有一个矢量化解决方案，就像这样 -

np.sort(np.random.rand(10000,365).argpartition(50,axis=1)[:,:50])

让我们验证输出数据的一致性

In [209]: out = np.sort(np.random.rand(10000,365).argpartition(50,axis=1)[:,:50])

In [210]: count = np.bincount(out.ravel(), minlength=365)

In [211]: print count.min(), count.max()
1277 1466

看起来很统一！让我们进入主要业务，即获取性能数据。

运行时测试

方法-

# Original approach
def org_app():
    return np.array([np.sort(np.random.choice(365, 50, replace=False)) for _ in range(10000)])

# @Nils Werner's soln
def sort_random_choice():
    return np.sort([np.random.choice(365, 50, replace=False) for _ in range(10000)], axis=1)

# @Miriam Farber's soln
def random_permute():
    l = np.array([True]*50 + [False]*315)
    total = np.arange(1,366)
    return np.array([total[np.random.permutation(l)] for _ in range(10000)])

# Proposed in this post
def argpartition_sort(nrows=10000, maxc=365, ncols=50):
    return np.sort(np.random.rand(nrows,maxc).argpartition(ncols,axis=1)[:,:ncols])

# @ayhan's soln
def while_loop():
    arr = np.random.choice(365, 10000)
    for i in range(49):
        arr2 = np.random.choice(365, 10000)
        comp = (arr2 == arr)
        while comp.any():
            duplicate = comp if i==0 else comp.any(axis=0)
            arr2[duplicate] = np.random.choice(365, duplicate.sum())
            comp = (arr2 == arr)
        arr = np.vstack([arr, arr2])
    arr = arr.T
    arr.sort(axis=1)
    return arr

时间安排 -

In [44]: %timeit org_app()
    ...: %timeit sort_random_choice()
    ...: %timeit random_permute()
    ...: %timeit argpartition_sort()
    ...: %timeit while_loop()
    ...: 
1 loops, best of 3: 258 ms per loop
1 loops, best of 3: 232 ms per loop
10 loops, best of 3: 166 ms per loop
10 loops, best of 3: 79.9 ms per loop
10 loops, best of 3: 58.6 ms per loop

【讨论】：

【解决方案4】：

一种可能的优化是通过将sort 放在循环之外来向量化它：

sample_dates = np.sort([np.random.choice(365, 50, replace=False) for _ in range(10000)], axis=1)

【讨论】：