计算矩阵不使用嵌套的if和for循环?

Question

我试图制作一个矩阵A(200,200)，其中A[i][j] = A[j][i]遵循概率p的伯努利分布，如果元素在同一索引簇中(即。

0-49、50-99、100-149、150-199)，否则则为q。 到目前为止，我尝试像下面这样直接实现它，

import numpy as np
from scipy.stats import bernoulli
def link_weight(p=0.05, q=0.04):
A = np.zeros((200, 200))
for i in range(len(A)):
for j in range(len(A)):
if i == j:
A[i][j] = 0
elif i < j:
if i < 50:
if j < 50:
A[i][j] = bernoulli.rvs(p)
A[j][i] = bernoulli.rvs(p)
else:
A[i][j] = bernoulli.rvs(q)
A[j][i] = bernoulli.rvs(q)
if 50 <= i <= 99:
if 50 <= j <= 99:
A[i][j] = bernoulli.rvs(p)
A[j][i] = bernoulli.rvs(p)
else:
A[i][j] = bernoulli.rvs(q)
A[j][i] = bernoulli.rvs(q)
if 100 <= i <= 149:
if 100 <= j <= 149:
A[i][j] = bernoulli.rvs(p)
A[j][i] = bernoulli.rvs(p)
else:
A[i][j] = bernoulli.rvs(q)
A[j][i] = bernoulli.rvs(q)
if 150 <= i <= 199:
if 150 <= j <= 199:
A[i][j] = bernoulli.rvs(p)
A[j][i] = bernoulli.rvs(p)
else:
A[i][j] = bernoulli.rvs(q)
A[j][i] = bernoulli.rvs(q)
return A

但是我想让它更容易读，更快。我如何通过避免嵌套if/for循环来证明这段代码?

Answer 1

Numpy索引非常强大，所以你可以使用它:

import numpy as np
from scipy.stats import bernoulli
def link_weight(p=0.05, q=0.04):
# Create bernouli array
A = bernoulli.rvs(q, size=(200, 200))
# Randomize values p=p
A[0:50, 0:50] = bernoulli.rvs(p, size=(50,50))
A[50:100, 50:100] = bernoulli.rvs(p, size=(50,50))
A[100:150, 100:150] = bernoulli.rvs(p, size=(50,50))
A[150:200, 150:200] = bernoulli.rvs(p, size=(50,50))
# Make symmetric
A = A + A.T - np.diag(A.diagonal())
return A
print(link_weight())

• 对称矩阵:https://stackoverflow.com/a/10806947/14536215
• 索引:https://numpy.org/doc/stable/reference/arrays.indexing.html

编辑:修改对称性使用:https://stackoverflow.com/a/2573982/14536215

Answer 2

为每组索引创建10个正方形数组:p为4,q为6。

a00 = bernoulli.rvs(p, size=(50, 50))
a11 = bernoulli.rvs(p, size=(50, 50))
a22 = bernoulli.rvs(p, size=(50, 50))
a33 = bernoulli.rvs(p, size=(50, 50))
a01 = bernoulli.rvs(q, size=(50, 50))
a02 = bernoulli.rvs(q, size=(50, 50))
a03 = bernoulli.rvs(q, size=(50, 50))
a12 = bernoulli.rvs(q, size=(50, 50))
a13 = bernoulli.rvs(q, size=(50, 50))
a23 = bernoulli.rvs(q, size=(50, 50))

(当然，这可以自动完成)

现在，将这些块组合成一个大矩阵:

A = np.vstack([np.hstack([a00, a01, a02, a03]),
np.hstack([a01, a11, a12, a13]),
np.hstack([a02, a12, a22, a23]),
np.hstack([a03, a13, a23, a33])])

最后，通过结合下三角形部分和它的转置，减去主对角线，使该数组对称:

np.tril(A) + np.tril(A, k=-1).T