Подтвердить что ты не робот

Как эффективно рассчитать стандартное отклонение?

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

[0] (0.01, 0.01, 0.02, 0.04, 0.03)
[1] (0.00, 0.02, 0.02, 0.03, 0.02)
[2] (0.01, 0.02, 0.02, 0.03, 0.02)
     ...
[n] (0.01, 0.00, 0.01, 0.05, 0.03)

Я хотел бы эффективно рассчитать среднее и стандартное отклонение для каждого индекса списка по всем элементам массива.

Чтобы сделать среднее значение, я перебрал массив и суммировал значение по заданному индексу списка. В конце я делю каждое значение в моем "списке средних" на n (я работаю с населением, а не с выборкой из населения).

Чтобы сделать стандартное отклонение, я повторяю цикл снова, теперь, когда у меня есть вычисленное среднее значение.

Я хотел бы избежать прохождения массива дважды, один раз для среднего значения, а затем один раз для SD (после того, как у меня есть среднее значение).

Существует ли эффективный метод для вычисления обоих значений, проходящий через массив только один раз? Подойдет любой код на интерпретируемом языке (например, Perl или Python) или псевдокод.

4b9b3361

ОТВЕТЫ

Ответ 1

Ответ заключается в использовании алгоритма Уэлфорда, который очень четко определен после "наивных методов" в:

Он более численно стабилен, чем двухпроходные или онлайновые простые сборщики квадратов, предложенные в других ответах. Стабильность действительно имеет значение только тогда, когда у вас есть много значений, близких друг к другу, поскольку они приводят к так называемому " катастрофическому аннулированию " в литературе с плавающей запятой.

Возможно, вы также захотите освежить разницу между делением на количество выборок (N) и N-1 в расчете дисперсии (квадратическое отклонение). Деление на N-1 приводит к непредвзятой оценке дисперсии от выборки, тогда как деление на N в среднем недооценивает дисперсию (поскольку она не учитывает дисперсию между средним по выборке и истинным средним).

Я написал две записи в блоге на эту тему, в которых более подробно рассказывается, как удалить предыдущие значения в Интернете:

Вы также можете взглянуть на мою Java-реализацию; javadoc, source и unit-тесты все онлайн:

Ответ 2

Основным ответом является накопление суммы как x (назовем его sum_x1), так и x 2 (назовите его 'sum_x2') по ходу. Значение стандартного отклонения:

stdev = sqrt((sum_x2 / n) - (mean * mean))

где

mean = sum_x / n

Это стандартное отклонение выборки; вы получаете стандартное отклонение населения, используя "n" вместо "n - 1" в качестве делителя.

Возможно, вам придется беспокоиться о численной стабильности при разнице между двумя большими числами, если вы имеете дело с большими образцами. Перейдите к внешним ссылкам в других ответах (Wikipedia и т.д.) Для получения дополнительной информации.

Ответ 3

Возможно, не то, что вы просили, но... Если вы используете массив numpy, он будет эффективно работать для вас:

from numpy import array

nums = array(((0.01, 0.01, 0.02, 0.04, 0.03),
              (0.00, 0.02, 0.02, 0.03, 0.02),
              (0.01, 0.02, 0.02, 0.03, 0.02),
              (0.01, 0.00, 0.01, 0.05, 0.03)))

print nums.std(axis=1)
# [ 0.0116619   0.00979796  0.00632456  0.01788854]

print nums.mean(axis=1)
# [ 0.022  0.018  0.02   0.02 ]

Кстати, в этом блоге есть интересная дискуссия и комментарии к однопроходным методам вычисления средств и отклонений:

Ответ 4

Вот буквальный чистый перевод Python реализации алгоритма Welford из http://www.johndcook.com/standard_deviation.html:

https://github.com/liyanage/python-modules/blob/master/running_stats.py

class RunningStats:

    def __init__(self):
        self.n = 0
        self.old_m = 0
        self.new_m = 0
        self.old_s = 0
        self.new_s = 0

    def clear(self):
        self.n = 0

    def push(self, x):
        self.n += 1

        if self.n == 1:
            self.old_m = self.new_m = x
            self.old_s = 0
        else:
            self.new_m = self.old_m + (x - self.old_m) / self.n
            self.new_s = self.old_s + (x - self.old_m) * (x - self.new_m)

            self.old_m = self.new_m
            self.old_s = self.new_s

    def mean(self):
        return self.new_m if self.n else 0.0

    def variance(self):
        return self.new_s / (self.n - 1) if self.n > 1 else 0.0

    def standard_deviation(self):
        return math.sqrt(self.variance())

Использование:

rs = RunningStats()
rs.push(17.0);
rs.push(19.0);
rs.push(24.0);

mean = rs.mean();
variance = rs.variance();
stdev = rs.standard_deviation();

Ответ 5

Модуль Python runstats Module предназначен именно для этого. Установить runstats из PyPI:

pip install runstats

Резюме runstats могут вызывать среднее, дисперсию, стандартное отклонение, асимметрию и эксцесс за один проход данных. Мы можем использовать это для создания вашей "бегущей" версии.

from runstats import Statistics

stats = [Statistics() for num in range(len(data[0]))]

for row in data:

    for index, val in enumerate(row):
        stats[index].push(val)

    for index, stat in enumerate(stats):
        print 'Index', index, 'mean:', stat.mean()
        print 'Index', index, 'standard deviation:', stat.stddev()

Статистические сводки основаны на методе Кнута и Велфорда для вычисления стандартного отклонения за один проход, как описано в Art of Computer Programming, Vol. 2, p. 232, 3-е издание. Преимущество этого - это численно стабильные и точные результаты.

Отказ от ответственности: Я являюсь автором модуля runstats Python.

Ответ 6

Statistics::Descriptive - очень приличный модуль Perl для этих типов вычислений:

#!/usr/bin/perl

use strict; use warnings;

use Statistics::Descriptive qw( :all );

my $data = [
    [ 0.01, 0.01, 0.02, 0.04, 0.03 ],
    [ 0.00, 0.02, 0.02, 0.03, 0.02 ],
    [ 0.01, 0.02, 0.02, 0.03, 0.02 ],
    [ 0.01, 0.00, 0.01, 0.05, 0.03 ],
];

my $stat = Statistics::Descriptive::Full->new;
# You also have the option of using sparse data structures

for my $ref ( @$data ) {
    $stat->add_data( @$ref );
    printf "Running mean: %f\n", $stat->mean;
    printf "Running stdev: %f\n", $stat->standard_deviation;
}
__END__

Вывод:

C:\Temp> g
Running mean: 0.022000
Running stdev: 0.013038
Running mean: 0.020000
Running stdev: 0.011547
Running mean: 0.020000
Running stdev: 0.010000
Running mean: 0.020000
Running stdev: 0.012566

Ответ 7

Посмотрите PDL (произносится как "piddle!" ).

Это язык данных Perl, который предназначен для высокоточной математики и научных вычислений.

Вот пример использования ваших цифр....

use strict;
use warnings;
use PDL;

my $figs = pdl [
    [0.01, 0.01, 0.02, 0.04, 0.03],
    [0.00, 0.02, 0.02, 0.03, 0.02],
    [0.01, 0.02, 0.02, 0.03, 0.02],
    [0.01, 0.00, 0.01, 0.05, 0.03],
];

my ( $mean, $prms, $median, $min, $max, $adev, $rms ) = statsover( $figs );

say "Mean scores:     ", $mean;
say "Std dev? (adev): ", $adev;
say "Std dev? (prms): ", $prms;
say "Std dev? (rms):  ", $rms;


Что производит:

Mean scores:     [0.022 0.018 0.02 0.02]
Std dev? (adev): [0.0104 0.0072 0.004 0.016]
Std dev? (prms): [0.013038405 0.010954451 0.0070710678 0.02]
Std dev? (rms):  [0.011661904 0.009797959 0.0063245553 0.017888544]


Посмотрите PDL:: Primitive для получения дополнительной информации о функции statsover. Это, по-видимому, свидетельствует о том, что ADEV является "стандартным отклонением".

Однако это может быть PRMS (что показывает Sinan Statistics:: Descriptive example) или RMS (что показывает пример NumPy). Я думаю, что один из этих трех должен быть прав: -)

Для получения дополнительной информации о PDL посмотрите:

Ответ 8

Насколько велик ваш массив? Если это не будет длинные элементы, не беспокойтесь о том, чтобы прокручивать его дважды. Код прост и легко протестирован.

Моим преимуществом было бы использовать расширение numpy, чтобы преобразовать массив массивов в двумерный массив numpy и получить стандартное отклонение:

>>> x = [ [ 1, 2, 4, 3, 4, 5 ], [ 3, 4, 5, 6, 7, 8 ] ] * 10
>>> import numpy
>>> a = numpy.array(x)
>>> a.std(axis=0) 
array([ 1. ,  1. ,  0.5,  1.5,  1.5,  1.5])
>>> a.mean(axis=0)
array([ 2. ,  3. ,  4.5,  4.5,  5.5,  6.5])

Если это не вариант и вам нужно чистое решение Python, продолжайте читать...

Если ваш массив

x = [ 
      [ 1, 2, 4, 3, 4, 5 ],
      [ 3, 4, 5, 6, 7, 8 ],
      ....
]

Тогда стандартное отклонение:

d = len(x[0])
n = len(x)
sum_x = [ sum(v[i] for v in x) for i in range(d) ]
sum_x2 = [ sum(v[i]**2 for v in x) for i in range(d) ]
std_dev = [ sqrt((sx2 - sx**2)/N)  for sx, sx2 in zip(sum_x, sum_x2) ]

Если вы настроены на циклическое перемещение массива только один раз, текущие суммы могут быть объединены.

sum_x  = [ 0 ] * d
sum_x2 = [ 0 ] * d
for v in x:
   for i, t in enumerate(v):
   sum_x[i] += t
   sum_x2[i] += t**2

Это не так элегантно, как решение для понимания списка выше.

Ответ 11

n=int(raw_input("Enter no. of terms:"))

L=[]

for i in range (1,n+1):

    x=float(raw_input("Enter term:"))

    L.append(x)

sum=0

for i in range(n):

    sum=sum+L[i]

avg=sum/n

sumdev=0

for j in range(n):

    sumdev=sumdev+(L[j]-avg)**2

dev=(sumdev/n)**0.5

print "Standard deviation is", dev

Ответ 12

Как говорится в следующем ответе: Предоставляет ли pandas/scipy/numpy функцию кумулятивного стандартного отклонения? Модуль Python Pandas содержит метод вычисления текущего или совокупного стандартного отклонения. Для этого вам нужно будет преобразовать ваши данные в фреймворк Pandas (или серию, если это 1D), но для этого есть функции.

Ответ 13

Здесь "однострочный", распространяемый по нескольким строкам, в стиле функционального программирования:

def variance(data, opt=0):
    return (lambda (m2, i, _): m2 / (opt + i - 1))(
        reduce(
            lambda (m2, i, avg), x:
            (
                m2 + (x - avg) ** 2 * i / (i + 1),
                i + 1,
                avg + (x - avg) / (i + 1)
            ),
            data,
            (0, 0, 0)))

Ответ 14

Мне нравится выражать обновление следующим образом:

def running_update(x, N, mu, var):
    '''
        @arg x: the current data sample
        @arg N : the number of previous samples
        @arg mu: the mean of the previous samples
        @arg var : the variance over the previous samples
        @retval (N+1, mu', var') -- updated mean, variance and count
    '''
    N = N + 1
    rho = 1.0/N
    d = x - mu
    mu += rho*d
    var += rho*((1-rho)*d**2 - var)
    return (N, mu, var)

так что однопроходная функция будет выглядеть так:

def one_pass(data):
    N = 0
    mu = 0.0
    var = 0.0
    for x in data:
        N = N + 1
        rho = 1.0/N
        d = x - mu
        mu += rho*d
        var += rho*((1-rho)*d**2 - var)
        # could yield here if you want partial results
   return (N, mu, var)

обратите внимание, что это вычисление выборочной дисперсии (1/N), а не объективная оценка дисперсии населения (которая использует коэффициент нормализации 1/(N-1)). В отличие от других ответов, переменная var, которая отслеживает текущую дисперсию, не увеличивается пропорционально количеству выборок. Во все времена это просто дисперсия множества выборок, которые мы видели до сих пор (нет окончательного "деления на n" при получении дисперсии).

В классе это будет выглядеть так:

class RunningMeanVar(object):
    def __init__(self):
        self.N = 0
        self.mu = 0.0
        self.var = 0.0
    def push(self, x):
        self.N = self.N + 1
        rho = 1.0/N
        d = x-self.mu
        self.mu += rho*d
        self.var += + rho*((1-rho)*d**2-self.var)
    # reset, accessors etc. can be setup as you see fit

Это также работает для взвешенных образцов:

def running_update(w, x, N, mu, var):
    '''
        @arg w: the weight of the current sample
        @arg x: the current data sample
        @arg mu: the mean of the previous N sample
        @arg var : the variance over the previous N samples
        @arg N : the number of previous samples
        @retval (N+w, mu', var') -- updated mean, variance and count
    '''
    N = N + w
    rho = w/N
    d = x - mu
    mu += rho*d
    var += rho*((1-rho)*d**2 - var)
    return (N, mu, var)