我正在尝试找到一种性能良好的方法来计算沿Numpy数组的轴距质心/重心的标准偏差。
在公式中,这是(很抱歉未对准):
我能想到的最好的办法是:
def weighted_com(A, axis, weights):
average = np.average(A, axis=axis, weights=weights)
return average * weights.sum() / A.sum(axis=axis).astype(float)
def weighted_std(A, axis):
weights = np.arange(A.shape[axis])
w1com2 = weighted_com(A, axis, weights)**2
w2com1 = weighted_com(A, axis, weights**2)
return np.sqrt(w2com1 - w1com2)
在
weighted_com
中,我需要校正权重和值的归一化(我猜这是一个丑陋的解决方法)。 weighted_std
可能很好。为了避免XY问题,我仍然要求我真正想要的(更好的
weighted_std
)而不是我的weighted_com
的更好版本。.astype(float)
是一种安全措施,因为我会将其应用于包含整数的直方图,当不在Python 3中或from __future__ import division
未激活时,由于整数除法会导致问题。 最佳答案
您想获取向量[1, 2, 3, ..., n]
的均值,方差和标准偏差-其中n
是输入矩阵A
沿感兴趣轴的维数-权重由矩阵A
本身给出。
具体来说,假设您要考虑沿垂直轴(axis=0
)的这些质量中心统计信息-这就是您编写的公式所对应的内容。对于固定列j
,您可以
n = A.shape[0]
r = np.arange(1, n+1)
mu = np.average(r, weights=A[:,j])
var = np.average(r**2, weights=A[:,j]) - mu**2
std = np.sqrt(var)
为了将不同列的所有计算放在一起,您必须堆叠一堆
r
副本(每列一个)以形成一个矩阵(我在下面的代码中称为R
)。稍加注意,您就可以使axis=0
和axis=1
都起作用。import numpy as np
def com_stats(A, axis=0):
A = A.astype(float) # if you are worried about int vs. float
n = A.shape[axis]
m = A.shape[(axis-1)%2]
r = np.arange(1, n+1)
R = np.vstack([r] * m)
if axis == 0:
R = R.T
mu = np.average(R, axis=axis, weights=A)
var = np.average(R**2, axis=axis, weights=A) - mu**2
std = np.sqrt(var)
return mu, var, std
例如,
A = np.array([[1, 1, 0], [1, 2, 1], [1, 1, 1]])
print(A)
# [[1 1 0]
# [1 2 1]
# [1 1 1]]
print(com_stats(A))
# (array([ 2. , 2. , 2.5]), # centre-of-mass mean by column
# array([ 0.66666667, 0.5 , 0.25 ]), # centre-of-mass variance by column
# array([ 0.81649658, 0.70710678, 0.5 ])) # centre-of-mass std by column
编辑:
可以避免使用
r
创建R
的内存副本来构建numpy.lib.stride_tricks
:交换行R = np.vstack([r] * m)
以上
from numpy.lib.stride_tricks import as_strided
R = as_strided(r, strides=(0, r.itemsize), shape=(m, n))
生成的
R
是一个(条纹的)ndarray
,其基础数组与r
的相同-绝对不会复制任何值。from numpy.lib.stride_tricks import as_strided
FMT = '''\
Shape: {}
Strides: {}
Position in memory: {}
Size in memory (bytes): {}
'''
def find_base_nbytes(obj):
if obj.base is not None:
return find_base_nbytes(obj.base)
return obj.nbytes
def stats(obj):
return FMT.format(obj.shape,
obj.strides,
obj.__array_interface__['data'][0],
find_base_nbytes(obj))
n=10
m=1000
r = np.arange(1, n+1)
R = np.vstack([r] * m)
S = as_strided(r, strides=(0, r.itemsize), shape=(m, n))
print(stats(r))
print(stats(R))
print(stats(S))
输出:
Shape: (10,)
Strides: (8,)
Position in memory: 4299744576
Size in memory (bytes): 80
Shape: (1000, 10)
Strides: (80, 8)
Position in memory: 4304464384
Size in memory (bytes): 80000
Shape: (1000, 10)
Strides: (0, 8)
Position in memory: 4299744576
Size in memory (bytes): 80
记下this SO answer和this one以获得有关如何获取跨越的
ndarray
的基础数组的内存地址和大小的说明。关于python - 沿Numpy阵列轴的质心的标准偏差,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/38556297/