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

Question

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

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

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

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

正如 cmets 中所指出的，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）。但是，由于您总是希望从三个元素中获取两个元素，因此您可以对要删除的元素进行采样（使用倒置概率），然后保留另外两个。这个 sn-p 做了这样的事情：

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)，但不是 @987654329 @、(2, 0) 或 (2, 1))。

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