有一种更快的方法来重复一段代码x次并取平均吗？

Question

开始于：

    a,b=np.ogrid[0:n+1:1,0:n+1:1]
    B=np.exp(1j*(np.pi/3)*np.abs(a-b))
    B[z,b] = np.exp(1j * (np.pi/3) * np.abs(z - b +x))
    B[a,z] = np.exp(1j * (np.pi/3) * np.abs(a - z +x))
    B[diag,diag]=1-1j/np.sqrt(3)

这将生成一个用作矩阵的n*n网格。你知道吗

n只是一个用来表示指数的数字，即a*b矩阵，其中a和b都上升到n

其中z是常数，我选择用B[z，B]和B[a，z]公式替换行和列。（基本上是相同的公式，但在np.abs公司（a-b））

矩阵的对角线由底线给出：

B[diag,diag]=1-1j/np.sqrt(3)

在哪里

diag=np.arange(n+1)

我想重复这段代码50次，其中唯一的变化是x，所以我将以50个版本的B结束np.ogrid公司. x是一个随机生成的数字，每次介于-0.8和0.8之间。你知道吗

x=np.random.uniform(-0.8,0.8)

我想生成50个版本的B，每次随机值为x，并使用以下定义取50个版本B的几何平均值：

def geo_mean(y):
   y = np.asarray(y)
   return np.prod(y ** (1.0 / y.shape[0]), axis=-1)

我尝试将B设置为某个索引的函数，然后使用a for uin range（）：循环，但这不起作用。除了复制和粘贴块50次，并表示每一个为B1，B2，B3等，我想不出另一种方法来解决这个问题。你知道吗

编辑：

我现在正在使用给定解决方案的一部分，以便清楚地显示我要查找的内容：

#A matrix with 50 random values between -0.8 and 0.8 to be used in the loop
X=np.random.uniform(-0.8,0.8, (50,1))
#constructing the base array before modification by random x values in position z
 a,b = np.ogrid[0:n+1:1,0:n+1:1]
 B = np.exp(1j * ( np.pi / 3) * np.abs( a - b ))
 B[diag,diag] = 1 - 1j / np.sqrt(3)
#list to store all modified arrays
 randomarrays = []
 for i in range( 0,50 ):
#copy array and modify it
    Bnew = np.copy( B )
    Bnew[z, b] = np.exp( 1j * ( np.pi / 3 ) * np.abs(z - b + X[i]))
    Bnew[a, z] = np.exp( 1j * ( np.pi / 3 ) * np.abs(a - z + X[i]))
    randomarrays.append(Bnew)

 Bstack = np.dstack(randomarrays)
#calculate the geometric mean value along the axis that was the row in 2D arrays
 B0 = geo_mean(Bstack)

在这个例子中，i的每一次迭代都使用相同的X值，我似乎找不到一种方法让I的每个新循环使用矩阵X中的下一个值。我不确定python中的++操作，我知道它在python中不起作用，我只是不知道如何使用python等价物。我希望一个循环使用一个X的值，然后下一个循环使用下一个值，依此类推，这样我就可以dstack所有矩阵的末尾，并为堆叠矩阵中的每个元素找到一个geo_mean。你知道吗

Answer 1

一种可行的方法是使用列表理解或生成器表达式：

>>> def f(n, z, x):
...     diag = np.arange(n+1)
...     a,b=np.ogrid[0:n+1:1,0:n+1:1]
...     B=np.exp(1j*(np.pi/3)*np.abs(a-b))
...     B[z,b] = np.exp(1j * (np.pi/3) * np.abs(z - b +x))
...     B[a,z] = np.exp(1j * (np.pi/3) * np.abs(a - z +x))
...     B[diag,diag]=1-1j/np.sqrt(3)
...     return B
... 
>>> X = np.random.uniform(-0.8, 0.8, (10,))
>>> np.prod((*map(np.power, map(f, 10*(4,), 10*(2,), X), 10 * (1/10,)),), axis=0)

但在你的具体例子中，我们可以做得更好； 利用恒等式exp(a) x exp(b) = exp(a + b)，我们可以将指数化后的几何平均数转换为指数化前的算术平均数。由于复数n次方根的多值性，几何平均值中出现了复数n次方根，因此需要小心一点。在下面的代码中，我们规范化了出现在range-pi，pi上的角度，以便始终与第n根相同。你知道吗

另外请注意，您提供的geo_mean函数肯定是错误的。它没有通过基本的健全性检查，即对同一事物的副本取平均值应该返回同一事物。我提供了一个更好的版本。它仍然不是完美的，但我认为实际上没有完美的解决方案，因为复根的非均匀性。你知道吗

因此，我建议在求幂之前取平均值。只要你的随机分布小于π，这就允许一个定义良好的平均过程，平均值实际上接近样本

import numpy as np

def f(n, z, X, do_it_pps_way=True):
    X = np.asanyarray(X)
    diag = np.arange(n+1)
    a,b=np.ogrid[0:n+1:1,0:n+1:1]
    B=np.exp(1j*(np.pi/3)*np.abs(a-b))
    X = X.reshape(-1,1,1)
    if do_it_pps_way:
        zbx = np.mean(np.abs(z-b+X), axis=0)
        azx = np.mean(np.abs(a-z+X), axis=0)
    else:
        zbx = np.mean((np.abs(z-b+X)+3) % 6 - 3, axis=0)
        azx = np.mean((np.abs(a-z+X)+3) % 6 - 3, axis=0)
    B[z,b] = np.exp(1j * (np.pi/3) * zbx)
    B[a,z] = np.exp(1j * (np.pi/3) * azx)
    B[diag,diag]=1-1j/np.sqrt(3)
    return B

def geo_mean(y):
    y = np.asarray(y)
    dim = len(y.shape)
    y = np.atleast_2d(y)
    v = np.prod(y, axis=0) ** (1.0 / y.shape[0])
    return v[0] if dim == 1 else v

def geo_mean_correct(y):
    y = np.asarray(y)
    return np.prod(y ** (1.0 / y.shape[0]), axis=0)

# demo that orig geo_mean is wrong
B = np.exp(1j * np.random.random((5, 5)))
# the mean of four times the same thing should be the same thing:
if not np.allclose(B, geo_mean([B, B, B, B])):
    print('geo_mean failed')
if np.allclose(B, geo_mean_correct([B, B, B, B])):
    print('but geo_mean_correct works')

n, z, m = 10, 3, 50

X = np.random.uniform(-0.8, 0.8, (m,))
B0 = f(n, z, X, do_it_pps_way=False)
B1 = np.prod((*map(np.power, map(f, m*(n,), m*(z,), X), m * (1/m,)),), axis=0)
B2 = geo_mean_correct([f(n, z, x) for x in X])

# This is the recommended way:
B_recommended = f(n, z, X, do_it_pps_way=True)

print()
print(np.allclose(B1, B0))
print(np.allclose(B2, B1))

Answer 2

我认为在处理问题时，您应该更多地依赖numpy功能。我自己不是一个numpy专家，所以肯定还有改进的空间：

from scipy.stats import gmean

n = 2
z = 1
a = np.arange(n + 1).reshape(1, n + 1)
#constructing the base array before modification by random x values in position z
B = np.exp(1j * (np.pi / 3) * np.abs(a - a.T))
B[a, a] = 1 - 1j / np.sqrt(3)
#list to store all modified arrays
random_arrays = []
for _ in range(50):
    #generate random x value
    x=np.random.uniform(-0.8, 0.8)
    #copy array and modify it
    B_new = np.copy(B)
    B_new[z, a] = np.exp(1j * (np.pi / 3) * np.abs(z - a + x))
    B_new[a, z] = np.exp(1j * (np.pi / 3) * np.abs(a - z + x))
    random_arrays.append(B_new)

#store all B arrays as a 3D array
B_stack = np.stack(random_arrays)
#calculate the geometric mean value along the axis that was the row in 2D arrays
geom_mean_for_rows = gmean(B_stack, axis = 2)

它使用来自^{}模块的几何平均函数来进行此计算。你知道吗