计算 numpy 数组之间的 MSE
Calculating MSE between numpy arrays
科学问题：
我有很多 3D 体积，它们中都有一个圆柱体，圆柱体 'upright' 在 z 轴上。包含圆柱体的体积非常嘈杂，就像超级嘈杂一样，你无法将圆柱体视为人类。如果我将这些体积中的 1000 个平均在一起，我可以看到圆柱体。每个卷都包含一个圆柱体的副本，但在某些情况下，圆柱体的方向可能不正确，所以我想要一种方法来解决这个问题。
我想出的解决方案：
我取了平均体积并将其投影到 z 轴和 x 轴（只是投影 numpy 数组），这样我在一个方向上得到一个漂亮的圆圈，在另一个方向上得到一个矩形。然后我获取每个 3D 体积并沿 Z 轴投射每个体积。 SNR 仍然很差，我看不到一个圆，但如果我对 2D 切片进行平均，我可以在平均几百个之后开始看到一个圆，并且在前 1000 个平均后很容易看到。为了计算每个体积的分数，我计算了 3D 体积的 MSE 向下投影 z 相对于其他三个阵列，第一个是向下投影 Z 的平均值，然后是向下投影 y 或 x 的平均值，最后是一个具有其中噪声的正态分布。
目前我有以下内容，其中 RawParticle 是 3D 数据，Ave 是平均值：
def normalise(array):
min = np.amin(array)
max = np.amax(array)
normarray = (array - min) / (max - min)
return normarray
def Noise(mag):
NoiseArray = np.random.normal(0, mag, size=(200,200,200))
return NoiseArray
#3D volume (normally use a for loop to iterate through al particles but for this example just showing one)
RawParticleProjected = np.sum(RawParticle, 0)
RawParticleProjectedNorm = normalise(RawParticleProjected)
#Average
AveProjected = np.sum(Ave, 0)
AveProjectedNorm = normalise(AveProjected)
#Noise Array
NoiseArray = Noise(0.5)
NoiseNorm = normalise(NoiseArray)
#Mean squared error
MSE = (np.square(np.subtract(RawParticleProjectedNorm, AveProjectedNorm))).mean()
然后我用平均轴 1 求和重复此操作，然后再次将原始粒子与噪声数组进行比较。
然而，当我比较应该都是圆的投影时，我的输出给出了最高的 MSE，如下所示：
我对 MSE 的理解是，其他两个群体应该具有高 MSE，而我同意的群体应该具有低 MSE。也许我的数据对于这种类型的分析来说太嘈杂了？但如果那是真的，那么我真的不知道如何做我正在做的事情。
如果有人能看一眼我的代码或启发我对 MSE 的理解，我将不胜感激。
感谢您花时间查看和阅读。
如果我正确理解了你的问题，你想计算出你的不同样本与平均值的接近程度。
通过比较样本，您希望找到包含迷失方向的圆柱体的异常值。
这非常符合 L2 norm 的定义，因此 MSE 应该在这里工作。
我会计算所有样本的平均 3D 图像，然后计算每个样本与该平均值的距离。然后我会比较这些值。
将样本与人工噪声图像进行比较的想法不错，但我不确定正态分布和标准化是否按计划进行。我可以是苹果和橘子。
我不认为沿着不同的轴观察投影是个好主意，
只是比较 3D 图像。
我用二维圆圈做了一些小测试，参数 alpha 表示图片中有多少噪音和多少圆圈。
（alpha=0 表示仅噪声，alpha=1 表示仅圆`）
import numpy as np
import matplotlib.pyplot as plt
grid_size = 20
radius = 5
mag = 1
def get_circle_stencil(radius):
xx, yy = np.meshgrid(np.linspace(-grid_size/2+1/2, grid_size/2-1/2, grid_size),
np.linspace(-grid_size/2+1/2, grid_size/2-1/2, grid_size))
dist = np.sqrt(xx**2 + yy**2)
inner = dist < (radius - 1/2)
return inner.astype(float)
def create_noise(mag, n_dim=2):
# return np.random.normal(0, mag, size=(grid_size,)*n_dim)
return np.random.uniform(0, mag, size=(grid_size,)*n_dim)
def create_noisy_sample(alpha, n_dim=2):
return (np.random.uniform(0, 1-alpha, size=(grid_size,)*n_dim) +
alpha*get_circle_stencil(radius))
fig = plt.figure()
ax = fig.subplots(nrows=3, ncols=3)
np.unravel_index(3, shape=(3, 3))
alpha_list = np.arange(9) / 10
for i, alpha in enumerate(alpha_list):
r, c = np.unravel_index(i, shape=(3, 3))
ax[r][c].imshow(*norm(create_noisy_sample(alpha=alpha)), cmap='Greys')
ax[r][c].set_title(f"alpha={alpha}")
ax[r][c].xaxis.set_ticklabels([])
ax[r][c].yaxis.set_ticklabels([])
然后我尝试了一些指标（mse、cosine similarity 和 binary cross entropy
并查看它们在不同 alpha 值下的表现。
def normalize(*args):
return [a / np.linalg.norm(a) for a in args]
def cosim(a, b):
return np.sum(a * b)
def mse(a, b):
return np.sqrt(np.sum((a-b)**2))
def bce(a, b):
# binary cross entropy implemented from tensorflow / keras
eps = 1e-7
res = a * np.log(b + eps)
res += (1 - a) * np.log(1 - b + eps)
return np.mean(-res)
我比较了NoiseA-NoiseB、Circle-Circle、Circle-Noise、Noise-Sample、Circle-Sample
alpha = 0.1
noise = create_noise(mag=1, grid_size=grid_size)
noise_b = create_noise(mag=1, grid_size=grid_size)
circle_reference = get_circle_stencil(radius=radius, grid_size=grid_size)
sample = create_noise(mag=1, grid_size=grid_size) + alpha * circle_reference
print('NoiseA-NoiseB:', mse(*norm(noise, noise_b))) # 0.718
print('Circle-Circle:', mse(*norm(circle, circle))) # 0.000
print('Circle-Noise:', mse(*norm(circle, noise))) # 1.168
print('Noise-Sample:', mse(*norm(noise, sample))) # 0.697
print('Circle-Sample:', mse(*norm(circle, sample))) # 1.100
print('NoiseA-NoiseB:', cosim(*norm(noise, noise_b))) # 0.741
print('Circle-Circle:', cosim(*norm(circle, circle))) # 1.000
print('Circle-Noise:', cosim(*norm(circle, noise))) # 0.317
print('Noise-Sample:', cosim(*norm(noise, sample))) # 0.757
print('Circle-Sample:', cosim(*norm(circle, sample))) # 0.393
print('NoiseA-NoiseB:', bce(*norm(noise, noise_b))) # 0.194
print('Circle-Circle:', bce(*norm(circle, circle))) # 0.057
print('Circle-Noise:', bce(*norm(circle, noise))) # 0.111
print('Noise-Circle:', bce(*norm(noise, circle))) # 0.636
print('Noise-Sample:', bce(*norm(noise, sample))) # 0.192
print('Circle-Sample:', bce(*norm(circle, sample))) # 0.104
n = 1000
ns = np.zeros(n)
cs = np.zeros(n)
for i, alpha in enumerate(np.linspace(0, 1, n)):
sample = create_noisy_sample(alpha=alpha)
ns[i] = mse(*norm(noise, sample))
cs[i] = mse(*norm(circle, sample))
fig, ax = plt.subplots()
ax.plot(np.linspace(0, 1, n), ns, c='b', label='noise-sample')
ax.plot(np.linspace(0, 1, n), cs, c='r', label='circle-sample')
ax.set_xlabel('alpha')
ax.set_ylabel('mse')
ax.legend()
对于您的问题，我只看比较 circle-sample（红色）。
不同的样本表现得好像它们具有不同的 alpha 值，您可以相应地对它们进行分组。而且您应该能够检测异常值，因为它们应该具有更高的 mse.
你说你必须组合 100-1000 张图片才能看到圆柱体，这表明你的问题中的 alpha 值非常小，但在 average mse
应该可以。
科学问题：
我有很多 3D 体积，它们中都有一个圆柱体，圆柱体 'upright' 在 z 轴上。包含圆柱体的体积非常嘈杂，就像超级嘈杂一样，你无法将圆柱体视为人类。如果我将这些体积中的 1000 个平均在一起，我可以看到圆柱体。每个卷都包含一个圆柱体的副本，但在某些情况下，圆柱体的方向可能不正确，所以我想要一种方法来解决这个问题。
我想出的解决方案：
我取了平均体积并将其投影到 z 轴和 x 轴（只是投影 numpy 数组），这样我在一个方向上得到一个漂亮的圆圈，在另一个方向上得到一个矩形。然后我获取每个 3D 体积并沿 Z 轴投射每个体积。 SNR 仍然很差，我看不到一个圆，但如果我对 2D 切片进行平均，我可以在平均几百个之后开始看到一个圆，并且在前 1000 个平均后很容易看到。为了计算每个体积的分数，我计算了 3D 体积的 MSE 向下投影 z 相对于其他三个阵列，第一个是向下投影 Z 的平均值，然后是向下投影 y 或 x 的平均值，最后是一个具有其中噪声的正态分布。
目前我有以下内容，其中 RawParticle 是 3D 数据，Ave 是平均值：
def normalise(array):
min = np.amin(array)
max = np.amax(array)
normarray = (array - min) / (max - min)
return normarray
def Noise(mag):
NoiseArray = np.random.normal(0, mag, size=(200,200,200))
return NoiseArray
#3D volume (normally use a for loop to iterate through al particles but for this example just showing one)
RawParticleProjected = np.sum(RawParticle, 0)
RawParticleProjectedNorm = normalise(RawParticleProjected)
#Average
AveProjected = np.sum(Ave, 0)
AveProjectedNorm = normalise(AveProjected)
#Noise Array
NoiseArray = Noise(0.5)
NoiseNorm = normalise(NoiseArray)
#Mean squared error
MSE = (np.square(np.subtract(RawParticleProjectedNorm, AveProjectedNorm))).mean()
然后我用平均轴 1 求和重复此操作，然后再次将原始粒子与噪声数组进行比较。
然而，当我比较应该都是圆的投影时，我的输出给出了最高的 MSE，如下所示：
我对 MSE 的理解是，其他两个群体应该具有高 MSE，而我同意的群体应该具有低 MSE。也许我的数据对于这种类型的分析来说太嘈杂了？但如果那是真的，那么我真的不知道如何做我正在做的事情。
如果有人能看一眼我的代码或启发我对 MSE 的理解，我将不胜感激。
感谢您花时间查看和阅读。
如果我正确理解了你的问题，你想计算出你的不同样本与平均值的接近程度。
通过比较样本，您希望找到包含迷失方向的圆柱体的异常值。
这非常符合 L2 norm 的定义，因此 MSE 应该在这里工作。
我会计算所有样本的平均 3D 图像，然后计算每个样本与该平均值的距离。然后我会比较这些值。
将样本与人工噪声图像进行比较的想法不错，但我不确定正态分布和标准化是否按计划进行。我可以是苹果和橘子。 我不认为沿着不同的轴观察投影是个好主意， 只是比较 3D 图像。
我用二维圆圈做了一些小测试，参数 alpha 表示图片中有多少噪音和多少圆圈。
（alpha=0 表示仅噪声，alpha=1 表示仅圆`）
import numpy as np
import matplotlib.pyplot as plt
grid_size = 20
radius = 5
mag = 1
def get_circle_stencil(radius):
xx, yy = np.meshgrid(np.linspace(-grid_size/2+1/2, grid_size/2-1/2, grid_size),
np.linspace(-grid_size/2+1/2, grid_size/2-1/2, grid_size))
dist = np.sqrt(xx**2 + yy**2)
inner = dist < (radius - 1/2)
return inner.astype(float)
def create_noise(mag, n_dim=2):
# return np.random.normal(0, mag, size=(grid_size,)*n_dim)
return np.random.uniform(0, mag, size=(grid_size,)*n_dim)
def create_noisy_sample(alpha, n_dim=2):
return (np.random.uniform(0, 1-alpha, size=(grid_size,)*n_dim) +
alpha*get_circle_stencil(radius))
fig = plt.figure()
ax = fig.subplots(nrows=3, ncols=3)
np.unravel_index(3, shape=(3, 3))
alpha_list = np.arange(9) / 10
for i, alpha in enumerate(alpha_list):
r, c = np.unravel_index(i, shape=(3, 3))
ax[r][c].imshow(*norm(create_noisy_sample(alpha=alpha)), cmap='Greys')
ax[r][c].set_title(f"alpha={alpha}")
ax[r][c].xaxis.set_ticklabels([])
ax[r][c].yaxis.set_ticklabels([])
然后我尝试了一些指标（mse、cosine similarity 和 binary cross entropy
并查看它们在不同 alpha 值下的表现。
def normalize(*args):
return [a / np.linalg.norm(a) for a in args]
def cosim(a, b):
return np.sum(a * b)
def mse(a, b):
return np.sqrt(np.sum((a-b)**2))
def bce(a, b):
# binary cross entropy implemented from tensorflow / keras
eps = 1e-7
res = a * np.log(b + eps)
res += (1 - a) * np.log(1 - b + eps)
return np.mean(-res)
我比较了NoiseA-NoiseB、Circle-Circle、Circle-Noise、Noise-Sample、Circle-Sample
alpha = 0.1
noise = create_noise(mag=1, grid_size=grid_size)
noise_b = create_noise(mag=1, grid_size=grid_size)
circle_reference = get_circle_stencil(radius=radius, grid_size=grid_size)
sample = create_noise(mag=1, grid_size=grid_size) + alpha * circle_reference
print('NoiseA-NoiseB:', mse(*norm(noise, noise_b))) # 0.718
print('Circle-Circle:', mse(*norm(circle, circle))) # 0.000
print('Circle-Noise:', mse(*norm(circle, noise))) # 1.168
print('Noise-Sample:', mse(*norm(noise, sample))) # 0.697
print('Circle-Sample:', mse(*norm(circle, sample))) # 1.100
print('NoiseA-NoiseB:', cosim(*norm(noise, noise_b))) # 0.741
print('Circle-Circle:', cosim(*norm(circle, circle))) # 1.000
print('Circle-Noise:', cosim(*norm(circle, noise))) # 0.317
print('Noise-Sample:', cosim(*norm(noise, sample))) # 0.757
print('Circle-Sample:', cosim(*norm(circle, sample))) # 0.393
print('NoiseA-NoiseB:', bce(*norm(noise, noise_b))) # 0.194
print('Circle-Circle:', bce(*norm(circle, circle))) # 0.057
print('Circle-Noise:', bce(*norm(circle, noise))) # 0.111
print('Noise-Circle:', bce(*norm(noise, circle))) # 0.636
print('Noise-Sample:', bce(*norm(noise, sample))) # 0.192
print('Circle-Sample:', bce(*norm(circle, sample))) # 0.104
n = 1000
ns = np.zeros(n)
cs = np.zeros(n)
for i, alpha in enumerate(np.linspace(0, 1, n)):
sample = create_noisy_sample(alpha=alpha)
ns[i] = mse(*norm(noise, sample))
cs[i] = mse(*norm(circle, sample))
fig, ax = plt.subplots()
ax.plot(np.linspace(0, 1, n), ns, c='b', label='noise-sample')
ax.plot(np.linspace(0, 1, n), cs, c='r', label='circle-sample')
ax.set_xlabel('alpha')
ax.set_ylabel('mse')
ax.legend()
对于您的问题，我只看比较 circle-sample（红色）。
不同的样本表现得好像它们具有不同的 alpha 值，您可以相应地对它们进行分组。而且您应该能够检测异常值，因为它们应该具有更高的 mse.
你说你必须组合 100-1000 张图片才能看到圆柱体，这表明你的问题中的 alpha 值非常小，但在 average mse
应该可以。