Стандартное отклонение от центра масс вдоль оси матрицы 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 чей базовый массив такой же, как rs - абсолютно никакого копирования каких-либо значений не происходит.

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,

Другие вопросы по тегам