python - 如何从随机加权矩阵中选择元素

Question

我对 python 还很陌生，并且在随机性方面遇到了一些问题。

我正在寻找类似的东西 RandomChoice在数学中。 我创建了一个维度为 10x3 的矩阵，随机数大于 0。让我们将每行的总和称为 s_i，其中 i=0,...,9

后来我想为每行选择 3 个元素中的 2 个(不重复)，加权概率为 s_ij/s_i 所以我需要这样的东西，但具有加权概率

n=10
aa=np.random.uniform(1000, 2500, (n,3))
print(aa)
help=[0,1,2]
dd=np.zeros((n,2))
for i in range(n):
    cc=random.sample(help,2)
    dd[i,0]=aa[i,cc[0]]
    dd[i,1]=aa[i,cc[1]]
print(dd)

在这里，另外速度是一个重要因素，因为我将在 Montecarlo 方法中使用它(这就是我从 Mathematica 切换到 Python 的原因)，我想，上面的代码可以得到很大的改进

提前感谢您的任何提示/帮助

编辑:我现在有以下内容，它可以工作，但对我来说看起来不太好

#pre-defined lists 
nn=3
aa=np.random.uniform(1000, 2500, (nn,3))
help1=[0,1,2]
help2=aa.sum(axis=1)
#now I create a weigthed prob list and fill it
help3=np.zeros((nn,3)) 
for i in range(nn):
    help3[i,0]=aa[i,0]/help2[i]
    help3[i,1]=aa[i,1]/help2[i]
    help3[i,2]=aa[i,2]/help2[i]
#every timestep when I have to choose 2 out of 3
help5=np.zeros((nn,2))
for i in range(nn):
    #cc=random.sample(help1,2)
    help4=np.random.choice(help1, 2, replace=False, p=[help3[i,0], help3[i,1], help3[i,2]])
    help5[i,0]=aa[i,cc[0]]
    help5[i,1]=aa[i,cc[1]]
print(help5)

Answer 1

正如评论中指出的，np.random.choice 接受权重参数，因此您可以简单地在循环中使用它:

import numpy as np

# Make input data
np.random.seed(0)
n = 10
aa = np.random.uniform(1000, 2500, (n, 3))
s = np.random.rand(n, 3)
# Normalize weights
s_norm = s / s.sum(1, keepdims=True)
# Output array
out = np.empty((n, 2), dtype=aa.dtype)
# Sample iteratively
for i in range(n):
    out[i] = aa[i, np.random.choice(3, size=2, replace=False, p=s_norm[i])]

但这并不是最有效的方法，因为通常使用向量化操作比循环快得多。不幸的是，我认为没有任何方法可以同时从多个分类分布中进行采样(请参阅 NumPy issue #15201 )。但是，由于您总是希望从三个元素中获取两个元素，因此您可以对要删除的元素进行采样(以倒置概率)，然后保留其他两个。这段代码做了类似的事情:

import numpy as np

# Make input data
np.random.seed(0)
n = 10
aa = np.random.uniform(1000, 2500, (n, 3))
s = np.random.rand(n, 3)
print(s)
# [[0.26455561 0.77423369 0.45615033]
#  [0.56843395 0.0187898  0.6176355 ]
#  [0.61209572 0.616934   0.94374808]
#  [0.6818203  0.3595079  0.43703195]
#  [0.6976312  0.06022547 0.66676672]
#  [0.67063787 0.21038256 0.1289263 ]
#  [0.31542835 0.36371077 0.57019677]
#  [0.43860151 0.98837384 0.10204481]
#  [0.20887676 0.16130952 0.65310833]
#  [0.2532916  0.46631077 0.24442559]]

# Invert weights
si = 1 / s
# Normalize
si_norm = si / si.sum(1, keepdims=True)
# Accumulate
si_cum = np.cumsum(si_norm, axis=1)
# Sample according to inverted probabilities
t = np.random.rand(n, 1)
idx = np.argmax(t < si_cum, axis=1)
# Get non-sampled indices
r = np.arange(3)
m = r != idx[:, np.newaxis]
choice = np.broadcast_to(r, m.shape)[m].reshape(n, -1)
print(choice)
# [[1 2]
#  [0 2]
#  [0 2]
#  [1 2]
#  [0 2]
#  [0 2]
#  [0 1]
#  [1 2]
#  [0 2]
#  [1 2]]

# Get corresponding data
out = np.take_along_axis(aa, choice, 1)

这样做的一个可能的缺点是所选元素将始终按顺序排列(也就是说，对于给定的行，您可能会得到索引对 (0, 1), (0, 2) 或 (1, 2)，但不是 (1, 0)、(2, 0)或 (2, 1))。

当然，如果您确实只需要一些样本，那么循环可能是最方便且可维护的解决方案，第二个只有在您需要大规模执行此操作时才有用。