Стандартное отклонение от центра масс вдоль оси матрицы Numpy
Я пытаюсь найти эффективный способ вычислить стандартное отклонение от центра масс / гравитации вдоль оси массива 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 ответ
Вы хотите взять среднее, дисперсию и стандартное отклонение вектора [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
s - абсолютно никакого копирования каких-либо значений не происходит.
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
Отдайте должное этому SO-ответу и этому объяснению о том, как получить адрес памяти и размер базового массива ndarray
,