Рассчитать коэффициент вариации окна в астропии - PullRequest
Новая
       173

Рассчитать коэффициент вариации окна в астропии

1 голос
/ 06 августа 2020

У меня есть массив, по которому я хочу рассчитать статистику использования астропии. У меня есть: 

from astropy.convolution import convolve
import numpy as np

x = np.random.randint(1, 10, size=(5, 5))
y = convolve(x, np.ones((3, 3)), boundary='extend', preserve_nan=True)

print(x)
print(y)

[[9 1 8 6 5]
 [4 2 1 8 4]
 [2 8 4 6 6]
 [8 4 8 5 6]
 [7 3 1 2 3]]

[[5.33333333 4.77777778 4.55555556 5.66666667 5.33333333]
 [4.55555556 4.33333333 4.88888889 5.33333333 5.55555556]
 [4.66666667 4.55555556 5.11111111 5.33333333 5.66666667]
 [5.44444444 5.         4.55555556 4.55555556 4.77777778]
 [6.         4.66666667 3.22222222 3.44444444 3.66666667]]

Каждый элемент в y - это среднее значение прямоугольника 3x3, нарисованного вокруг этой позиции в x. Вместо среднего я хотел бы иметь коэффициент вариации (стандартное отклонение, деленное на среднее значение). Я не уверен, можно ли это сделать в астропии или мне нужно использовать что-то еще, например 

from scipy.stats import variation

1 Ответ

2 голосов
/ 06 августа 2020

Astropy строится на numpy и scipy. Numpy - это библиотека массивов нижнего уровня, которая реализует хранение данных и базовые c операции, которые используются библиотеками более высокого уровня, такими как scipy и astropy. Понимание того, как работают массивы numpy, поможет вам работать с астропией.

Поскольку вы хотите вести статистику в скользящем окне, что-то вроде scipy.stats.variation не поможет вам напрямую. Он вычисляет статистику по элементам, сгруппированным по осям, а это не то, как сгруппированы ваши элементы. Вы можете смоделировать скользящее окно по осям, используя что-то вроде np.lib.stride_tricks.as_strided, но это опасно, подвержено ошибкам и не очень эффективно. Вместо этого я предлагаю использовать метод свертки, аналогичный тому, что вы уже сделали.

Во-первых, помните, что дисперсия может быть выражена как 

var = sum(x^2) / N - (sum(x) / N)^2

Итак, давайте сделаем это: 

kernel = np.ones((3, 3))
mean_x = convolve(x, kernel, boundary='extend', preserve_nan=True)
mean_square_x = convolve(x**2, kernel, boundary='extend', preserve_nan=True)

Теперь у вас есть sum(x) / N для каждого окна и sum(x^2) / N. Стандартное отклонение - это просто квадрат root дисперсии: 

std_x = np.sqrt(mean_square_x - mean_x**2)

Требуемый результат - это соотношение 

result = std_x / mean_x

Вы можете выполнить все вычисления более кратко, как 

m = convolve(x, kernel, boundary='extend', preserve_nan=True)
result = np.sqrt(convolve(x**2, kernel, boundary='extend', preserve_nan=True) - m**2) / m
...