Алгоритмы расчета дисперсии - Algorithms for calculating variance

Алгоритмы расчета дисперсии играть важную роль в вычислительная статистика. Ключевая трудность в дизайне хорошего алгоритмы для этой проблемы состоит в том, что формулы для отклонение может включать суммы квадратов, что может привести к числовая нестабильность а также арифметическое переполнение при работе с большими значениями.

Наивный алгоритм

Формула для расчета дисперсии всего численность населения размера N является:

${displaystyle sigma ^ {2} = {overline {(x ^ {2})}} - {ar {x}} ^ {2} = displaystyle {frac {sum _ {i = 1} ^ {N} x_ {i } ^ {2} - (сумма _ {i = 1} ^ {N} x_ {i}) ^ {2} / N} {N}}.!}$

С помощью Поправка Бесселя рассчитать беспристрастный оценка дисперсии совокупности от конечного образец из п наблюдений формула:

${displaystyle s ^ {2} = left ({frac {sum _ {i = 1} ^ {n} x_ {i} ^ {2}} {n}} - left ({frac {sum _ {i = 1}) ^ {n} x_ {i}} {n}} ight) ^ {2} ight) cdot {frac {n} {n-1}}.!}$

Следовательно, наивный алгоритм вычисления предполагаемой дисперсии имеет следующий вид:

Этот алгоритм можно легко адаптировать для вычисления дисперсии конечной совокупности: просто разделите на N вместо п - 1 в последней строке.

Потому что SumSq и (Сумма × Сумма) /п могут быть очень похожие числа, отмена может привести к точность результата будет намного меньше, чем присущая точность арифметика с плавающей запятой используется для выполнения вычислений. Таким образом, этот алгоритм не следует использовать на практике,^[1][2] и было предложено несколько альтернативных, численно устойчивых алгоритмов.^[3] Это особенно плохо, если стандартное отклонение мало по сравнению со средним значением. Однако алгоритм можно улучшить, приняв метод предполагаемое среднее.

Вычисление смещенных данных

Разница составляет инвариантный относительно изменений в параметр местоположения, свойство, которое можно использовать, чтобы избежать катастрофической отмены в этой формуле.

${displaystyle operatorname {Var} (X-K) = operatorname {Var} (X).}$

с ${displaystyle K}$ любая константа, которая приводит к новой формуле

${displaystyle s ^ {2} = displaystyle {frac {sum _ {i = 1} ^ {n} (x_ {i} -K) ^ {2} - (sum _ {i = 1} ^ {n} (x_ {i} -K)) ^ {2} / n} {n-1}}.!}$

чем ближе ${displaystyle K}$ соответствует среднему значению, тем точнее будет результат, но простой выбор значения в пределах диапазона выборки гарантирует желаемую стабильность. Если значения ${displaystyle (x_ {i} -K)}$ малы, то нет проблем с суммой его квадратов, наоборот, если они большие, это обязательно означает, что и дисперсия велика. В любом случае второй член в формуле всегда меньше первого, поэтому отмены не происходит.^[2]

Если в качестве ${displaystyle K}$ алгоритм можно записать на Язык программирования Python в качестве

Эта статья возможно содержит оригинальные исследования. Пожалуйста Улучши это к проверка заявленные претензии и добавление встроенные цитаты. Заявления, содержащие только оригинальные исследования, следует удалить. (Август 2019 г.) (Узнайте, как и когда удалить этот шаблон сообщения)

def shift_data_variance(данные):    если len(данные) < 2:        возвращаться 0.0    K = данные[0]    п = Бывший = Ex2 = 0.0    за Икс в данные:        п = п + 1        Бывший += Икс - K        Ex2 += (Икс - K) * (Икс - K)    отклонение = (Ex2 - (Бывший * Бывший) / п) / (п - 1)    # используйте n вместо (n-1), если хотите вычислить точную дисперсию заданных данных    # используйте (n-1), если данные являются выборками для большей совокупности    возвращаться отклонение

Эта формула также облегчает инкрементные вычисления, которые можно выразить как

K = п = Бывший = Ex2 = 0.0def add_variable(Икс):    Глобальный K, п, Бывший, Ex2    если п == 0:        K = Икс    п += 1    Бывший += Икс - K    Ex2 += (Икс - K) * (Икс - K)def remove_variable(Икс):    Глобальный K, п, Бывший, Ex2    п -= 1    Бывший -= Икс - K    Ex2 -= (Икс - K) * (Икс - K)def get_mean():    Глобальный K, п, Бывший    возвращаться K + Бывший / пdef get_variance():    Глобальный п, Бывший, Ex2    возвращаться (Ex2 - (Бывший * Бывший) / п) / (п - 1)

Двухпроходный алгоритм

Альтернативный подход, использующий другую формулу для дисперсии, сначала вычисляет выборочное среднее,

${displaystyle {ar {x}} = {frac {sum _ {j = 1} ^ {n} x_ {j}} {n}},}$

а затем вычисляет сумму квадратов отличий от среднего,

${displaystyle {ext {sample variance}} = s ^ {2} = displaystyle {frac {sum _ {i = 1} ^ {n} (x_ {i} - {ar {x}}) ^ {2}} { п-1}} ,!}$

куда s стандартное отклонение. Это дается следующим кодом:

def two_pass_variance(данные):    п = сумма1 = сумма2 = 0    за Икс в данные:        п += 1        сумма1 += Икс    иметь в виду = сумма1 / п    за Икс в данные:        сумма2 += (Икс - иметь в виду) * (Икс - иметь в виду)    отклонение = сумма2 / (п - 1)    возвращаться отклонение

Этот алгоритм численно устойчив, если п маленький.^[1][4] Однако результаты обоих этих простых алгоритмов («наивный» и «двухпроходный») могут чрезмерно зависеть от порядка данных и могут давать плохие результаты для очень больших наборов данных из-за повторяющейся ошибки округления при накоплении суммы. Такие методы, как компенсированное суммирование может быть использован для борьбы с этой ошибкой в ​​определенной степени.

Онлайн алгоритм Велфорда

Часто бывает полезно вычислить дисперсию за один проход, проверяя каждое значение. ${displaystyle x_ {i}}$ только однажды; например, когда данные собираются без достаточного объема памяти для хранения всех значений или когда затраты на доступ к памяти преобладают над затратами на вычисления. Для такого онлайн алгоритм, а отношение повторения требуется между величинами, из которых можно рассчитать требуемую статистику численно стабильным образом.

Следующие формулы можно использовать для обновления иметь в виду и (оценочная) дисперсия последовательности для дополнительного элемента Икс_п. Здесь, Икс_п обозначает выборочное среднее первого п образцы (Икс₁, ..., Икс_п), s²
_п их выборочная дисперсия, и σ²
_п их дисперсия населения.

${displaystyle {ar {x}} _ {n} = {frac {(n-1), {ar {x}} _ {n-1} + x_ {n}} {n}} = {ar {x} } _ {n-1} + {frac {x_ {n} - {ar {x}} _ {n-1}} {n}}!}$

${displaystyle s_ {n} ^ {2} = {frac {n-2} {n-1}}, s_ {n-1} ^ {2} + {frac {(x_ {n} - {ar {x}) } _ {n-1}) ^ {2}} {n}} = s_ {n-1} ^ {2} + {frac {(x_ {n} - {ar {x}} _ {n-1} ) ^ {2}} {n}} - {frac {s_ {n-1} ^ {2}} {n-1}}, quad n> 1}$

${displaystyle sigma _ {n} ^ {2} = {frac {(n-1), sigma _ {n-1} ^ {2} + (x_ {n} - {ar {x}} _ {n-1) }) (x_ {n} - {ar {x}} _ {n})} {n}} = sigma _ {n-1} ^ {2} + {frac {(x_ {n} - {ar {x }} _ {n-1}) (x_ {n} - {ar {x}} _ {n}) - sigma _ {n-1} ^ {2}} {n}}.}$

Эти формулы страдают числовой нестабильностью, поскольку они многократно вычитают небольшое число из большого числа, которое масштабируется с п. Лучшее количество для обновления - это сумма квадратов отличий от текущего среднего, ${displaystyle extstyle sum _ {i = 1} ^ {n} (x_ {i} - {ar {x}} _ {n}) ^ {2}}$, здесь обозначено ${displaystyle M_ {2, n}}$:

${displaystyle {egin {align} M_ {2, n} & = M_ {2, n-1} + (x_ {n} - {ar {x}} _ {n-1}) (x_ {n} - { ar {x}} _ {n}) [4pt] s_ {n} ^ {2} & = {frac {M_ {2, n}} {n-1}} [4pt] sigma _ {n} ^ {2} & = {frac {M_ {2, n}} {n}} конец {выровнено}}}$

Этот алгоритм был найден Велфордом,^[5][6] и это было тщательно проанализировано.^[2][7] Также принято обозначать ${displaystyle M_ {k} = {ar {x}} _ {k}}$ и ${displaystyle S_ {k} = M_ {2, k}}$.^[8]

Ниже приведен пример реализации алгоритма Велфорда на Python.

# Для нового значения newValue вычислить новый счетчик, новое среднее значение, новое значение M2.# mean накапливает среднее значение всего набора данных# M2 суммирует квадрат расстояния от среднего# count объединяет количество просмотренных образцовdef Обновить(existingAggregate, newValue):    (считать, иметь в виду, M2) = existingAggregate    считать += 1    дельта = newValue - иметь в виду    иметь в виду += дельта / считать    дельта2 = newValue - иметь в виду    M2 += дельта * дельта2    возвращаться (считать, иметь в виду, M2)# Получить среднее значение, дисперсию и дисперсию выборки из агрегатаdef завершить(existingAggregate):    (считать, иметь в виду, M2) = existingAggregate    если считать < 2:        возвращаться плавать("нан")    еще:        (иметь в виду, отклонение, выборочная дисперсия) = (иметь в виду, M2 / считать, M2 / (считать - 1))        возвращаться (иметь в виду, отклонение, выборочная дисперсия)

Этот алгоритм гораздо менее подвержен потере точности из-за катастрофическая отмена, но может быть не таким эффективным из-за операции деления внутри цикла. Для особенно надежного двухпроходного алгоритма вычисления дисперсии можно сначала вычислить и вычесть оценку среднего, а затем использовать этот алгоритм для остатков.

В параллельный алгоритм ниже показано, как объединить несколько наборов статистических данных, рассчитанных в режиме онлайн.

Взвешенный инкрементальный алгоритм

Алгоритм может быть расширен для обработки неравных весов выборок, заменив простой счетчик п с суммой видимых весов. Запад (1979)^[9] предлагает это инкрементальный алгоритм:

def weighted_incremental_variance(data_weight_pairs):    w_sum = w_sum2 = иметь в виду = S = 0    за Икс, ш в data_weight_pairs:  # В качестве альтернативы "для x, w в zip (данные, веса):"        w_sum = w_sum + ш        w_sum2 = w_sum2 + ш * ш        mean_old = иметь в виду        иметь в виду = mean_old + (ш / w_sum) * (Икс - mean_old)        S = S + ш * (Икс - mean_old) * (Икс - иметь в виду)    Population_variance = S / w_sum    # Поправка Бесселя для взвешенных выборок    # Частотные веса    sample_frequency_variance = S / (w_sum - 1)    # Веса надежности    sample_reliability_variance = S / (w_sum - w_sum2 / w_sum)

Параллельный алгоритм

Чан и др.^[10] обратите внимание, что онлайн-алгоритм Велфорда, описанный выше, является частным случаем алгоритма, который работает для объединения произвольных наборов ${displaystyle A}$ и ${displaystyle B}$:

${displaystyle {egin {align} n_ {AB} & = n_ {A} + n_ {B} delta & = {ar {x}} _ {B} - {ar {x}} _ {A} {ar {x}} _ {AB} & = {ar {x}} _ {A} + delta cdot {frac {n_ {B}} {n_ {AB}}} M_ {2, AB} & = M_ {2 , A} + M_ {2, B} + delta ^ {2} cdot {frac {n_ {A} n_ {B}} {n_ {AB}}} end {выровнено}}}$.

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

Метод Чана для оценки среднего численно неустойчив, когда ${displaystyle n_ {A} примерно n_ {B}}$ и оба они велики, потому что числовая ошибка в ${displaystyle delta = {ar {x}} _ {B} - {ar {x}} _ {A}}$ не уменьшается так, как в ${displaystyle n_ {B} = 1}$ дело. В таких случаях предпочитайте ${displaystyle {ar {x}} _ {AB} = {frac {n_ {A} {ar {x}} _ {A} + n_ {B} {ar {x}} _ {B}} {n_ {AB} }}}}$.

def parallel_variance(н_а, avg_a, M2_a, н_б, avg_b, M2_b):    п = н_а + н_б    дельта = avg_b - avg_a    M2 = M2_a + M2_b + дельта ** 2 * н_а * н_б / п    var_ab = M2 / (п - 1)    возвращаться var_ab

Это можно обобщить, чтобы разрешить распараллеливание с AVX, с GPU, и компьютерные кластеры, и ковариантности.^[3]

Пример

Предположим, что все операции с плавающей запятой используют стандартные IEEE 754 с двойной точностью арифметика. Рассмотрим выборку (4, 7, 13, 16) из бесконечной совокупности. На основе этой выборки оценочное среднее значение совокупности составляет 10, а несмещенная оценка дисперсии совокупности - 30. И наивный алгоритм, и двухпроходный алгоритм вычисляют эти значения правильно.

Далее рассмотрим образец (10⁸ + 4, 10⁸ + 7, 10⁸ + 13, 10⁸ + 16), что приводит к той же оценке дисперсии, что и в первом примере. Двухпроходный алгоритм правильно вычисляет эту оценку дисперсии, но наивный алгоритм возвращает 29,333333333333332 вместо 30.

Хотя эта потеря точности может быть терпимой и рассматриваться как незначительный недостаток наивного алгоритма, дальнейшее увеличение смещения делает ошибку катастрофической. Рассмотрим образец (10⁹ + 4, 10⁹ + 7, 10⁹ + 13, 10⁹ + 16). Опять же, оценочная дисперсия совокупности, равная 30, правильно вычисляется двухпроходным алгоритмом, но теперь простой алгоритм вычисляет ее как -170,66666666666666. Это серьезная проблема с наивным алгоритмом и связана с катастрофическая отмена при вычитании двух одинаковых чисел на заключительном этапе алгоритма.

Статистика высшего порядка

Terriberry^[11] расширяет формулы Чана на вычисление третьего и четвертого центральные моменты, например, при оценке перекос и эксцесс:

${displaystyle {egin {align} M_ {3, X} = M_ {3, A} + M_ {3, B} & {} + delta ^ {3} {frac {n_ {A} n_ {B} (n_ { A} -n_ {B})} {n_ {X} ^ {2}}} + 3delta {frac {n_ {A} M_ {2, B} -n_ {B} M_ {2, A}} {n_ { X}}} [6pt] M_ {4, X} = M_ {4, A} + M_ {4, B} & {} + delta ^ {4} {frac {n_ {A} n_ {B} left ( n_ {A} ^ {2} -n_ {A} n_ {B} + n_ {B} ^ {2} ight)} {n_ {X} ^ {3}}} [6pt] & {} + 6delta ^ {2} {frac {n_ {A} ^ {2} M_ {2, B} + n_ {B} ^ {2} M_ {2, A}} {n_ {X} ^ {2}}} + 4delta { гидроразрыв {n_ {A} M_ {3, B} -n_ {B} M_ {3, A}} {n_ {X}}} конец {выровнено}}}$

Здесь ${displaystyle M_ {k}}$ снова суммы степеней отличий от среднего ${displaystyle sum (x- {overline {x}}) ^ {k}}$, давая

${displaystyle {egin {align} & {ext {skewness}} = g_ {1} = {frac {{sqrt {n}} M_ {3}} {M_ {2} ^ {3/2}}}}, [ 4pt] & {ext {kurtosis}} = g_ {2} = {frac {nM_ {4}} {M_ {2} ^ {2}}} - 3.end {выровнено}}}$

Для инкрементального случая (т. Е. ${displaystyle B = {x}}$), это упрощает:

${displaystyle {egin {align} delta & = xm [5pt] m '& = m + {frac {delta} {n}} [5pt] M_ {2}' & = M_ {2} + delta ^ {2} {frac {n-1} {n}} [5pt] M_ {3} '& = M_ {3} + delta ^ {3} {frac {(n-1) (n-2)} {n ^ { 2}}} - {frac {3delta M_ {2}} {n}} [5pt] M_ {4} '& = M_ {4} + {frac {delta ^ {4} (n-1) (n ^ {2} -3n + 3)} {n ^ {3}}} + {frac {6delta ^ {2} M_ {2}} {n ^ {2}}} - {frac {4delta M_ {3}} { n}} конец {выровнен}}}$

Сохраняя ценность ${displaystyle delta / n}$, требуется только одна операция деления, и, таким образом, можно рассчитать статистику более высокого порядка с небольшими дополнительными затратами.

Пример онлайн-алгоритма эксцесса, реализованного, как описано, это:

def online_kurtosis(данные):    п = иметь в виду = M2 = M3 = M4 = 0    за Икс в данные:        n1 = п        п = п + 1        дельта = Икс - иметь в виду        delta_n = дельта / п        delta_n2 = delta_n * delta_n        термин1 = дельта * delta_n * n1        иметь в виду = иметь в виду + delta_n        M4 = M4 + термин1 * delta_n2 * (п*п - 3*п + 3) + 6 * delta_n2 * M2 - 4 * delta_n * M3        M3 = M3 + термин1 * delta_n * (п - 2) - 3 * delta_n * M2        M2 = M2 + термин1    # Обратите внимание: вы также можете рассчитать дисперсию, используя M2, и асимметрию, используя M3    эксцесс = (п * M4) / (M2 * M2) - 3    возвращаться эксцесс

Пеба́^[12]далее расширяет эти результаты до произвольного порядка центральные моменты, для инкрементального и попарного случаев, а впоследствии Пеба и др.^[13]для взвешенных и составных моментов. Там же можно найти аналогичные формулы для ковариация.

Чой и Свитмен^[14]предлагают два альтернативных метода вычисления асимметрии и эксцесса, каждый из которых может существенно сэкономить требования к памяти компьютера и процессорного времени в определенных приложениях. Первый подход заключается в вычислении статистических моментов путем разделения данных на ячейки и последующего вычисления моментов из геометрии полученной гистограммы, которая фактически становится однопроходный алгоритм для более высоких моментов. Одно из преимуществ состоит в том, что вычисления статистического момента могут выполняться с произвольной точностью, так что вычисления могут быть настроены с точностью, например, до формата хранения данных или исходного измерительного оборудования. Относительная гистограмма случайной величины может быть построена обычным способом: диапазон потенциальных значений делится на интервалы, а количество вхождений в каждом интервале подсчитывается и строится таким образом, чтобы площадь каждого прямоугольника равнялась части значений выборки в пределах эта корзина:

${displaystyle H (x_ {k}) = {гидроразрыв {h (x_ {k})} {A}}}$

куда ${displaystyle h (x_ {k})}$ и ${displaystyle H (x_ {k})}$ представляют частоту и относительную частоту в ячейке ${displaystyle x_ {k}}$ и ${displaystyle A = sum _ {k = 1} ^ {K} h (x_ {k}), Delta x_ {k}}$ - общая площадь гистограммы. После этой нормализации ${displaystyle n}$ сырые моменты и центральные моменты ${displaystyle x (t)}$ можно рассчитать по относительной гистограмме:

${displaystyle m_ {n} ^ {(h)} = sum _ {k = 1} ^ {K} x_ {k} ^ {n} H (x_ {k}), дельта x_ {k} = {frac {1 } {A}} сумма _ {k = 1} ^ {K} x_ {k} ^ {n} h (x_ {k}), Delta x_ {k}}$

${displaystyle heta _ {n} ^ {(h)} = sum _ {k = 1} ^ {K} {Big (} x_ {k} -m_ {1} ^ {(h)} {Big)} ^ { n}, H (x_ {k}), Delta x_ {k} = {frac {1} {A}} sum _ {k = 1} ^ {K} {Big (} x_ {k} -m_ {1} ^ {(h)} {Большой)} ^ {n} h (x_ {k}), Дельта x_ {k}}$

где верхний индекс ${displaystyle ^ {(h)}}$ указывает, что моменты рассчитываются по гистограмме. Для постоянной ширины бункера ${displaystyle Delta x_ {k} = Delta x}$ эти два выражения можно упростить, используя ${displaystyle I = A / Delta x}$:

${displaystyle m_ {n} ^ {(h)} = {frac {1} {I}} sum _ {k = 1} ^ {K} x_ {k} ^ {n}, h (x_ {k})}$

${displaystyle heta _ {n} ^ {(h)} = {frac {1} {I}} sum _ {k = 1} ^ {K} {Big (} x_ {k} -m_ {1} ^ {( h)} {Большой)} ^ {n} h (x_ {k})}$

Второй подход Чоя и Свитмена^[14] представляет собой аналитическую методологию для объединения статистических моментов из отдельных сегментов истории времени так, чтобы результирующие общие моменты были моментами полной истории времени. Эта методология может использоваться для параллельного вычисления статистических моментов с последующей комбинацией этих моментов или для комбинации статистических моментов, вычисляемых в последовательные моменты времени.

Если ${displaystyle Q}$ известны наборы статистических моментов:${displaystyle (gamma _ {0, q}, mu _ {q}, sigma _ {q} ^ {2}, alpha _ {3, q}, alpha _ {4, q}) quad}$ за ${displaystyle q = 1,2, ldots, Q}$, то каждый ${displaystyle gamma _ {n}}$ можно выразить через эквивалент ${displaystyle n}$ сырые моменты:

${displaystyle gamma _ {n, q} = m_ {n, q} gamma _ {0, q} qquad quad {extrm {for}} quad n = 1,2,3,4quad {ext {and}} quad q = 1,2, точки, Q}$

куда ${displaystyle gamma _ {0, q}}$ обычно принимается за продолжительность ${displaystyle q ^ {th}}$ история времени, или количество баллов, если ${displaystyle Delta t}$ постоянно.

Преимущество выражения статистических моментов в терминах ${displaystyle gamma}$ это то ${displaystyle Q}$ наборы могут быть объединены путем сложения, и нет верхнего предела значения ${displaystyle Q}$.

${displaystyle gamma _ {n, c} = sum _ {q = 1} ^ {Q} gamma _ {n, q} quad quad {ext {for}} n = 0,1,2,3,4}$

где нижний индекс ${displaystyle _ {c}}$ представляет собой объединенную временную историю или комбинированную ${displaystyle gamma}$. Эти комбинированные значения ${displaystyle gamma}$ затем могут быть обратно преобразованы в необработанные моменты, представляющие полную объединенную временную историю

${displaystyle m_ {n, c} = {frac {gamma _ {n, c}} {gamma _ {0, c}}} quad {ext {for}} n = 1,2,3,4}$

Известные отношения между необработанными моментами (${displaystyle m_ {n}}$) и центральные моменты (${displaystyle heta _ {n} = имя оператора {E} [(x-mu) ^ {n}])}$) затем используются для вычисления центральных моментов объединенной истории времени. Наконец, статистические моменты объединенной истории вычисляются из центральных моментов:

${displaystyle mu _ {c} = m_ {1, c} qquad sigma _ {c} ^ {2} = heta _ {2, c} qquad alpha _ {3, c} = {frac {heta _ {3, c }} {sigma _ {c} ^ {3}}} qquad alpha _ {4, c} = {frac {heta _ {4, c}} {sigma _ {c} ^ {4}}} - 3}$

Ковариация

Очень похожие алгоритмы можно использовать для вычисления ковариация.

Наивный алгоритм

Наивный алгоритм:

${displaystyle operatorname {Cov} (X, Y) = displaystyle {frac {sum _ {i = 1} ^ {n} x_ {i} y_ {i} - (sum _ {i = 1} ^ {n} x_ { i}) (sum _ {i = 1} ^ {n} y_ {i}) / n} {n}}.!}$

Для приведенного выше алгоритма можно использовать следующий код Python:

def наивная ковариация(data1, данные2):    п = len(data1)    сумма12 = 0    сумма1 = сумма(data1)    сумма2 = сумма(данные2)    за i1, i2 в застегивать(data1, данные2):        сумма12 += i1 * i2    ковариация = (сумма12 - сумма1 * сумма2 / п) / п    возвращаться ковариация

С оценкой среднего

Что касается дисперсии, ковариация двух случайных величин также инвариантна относительно сдвига, поэтому при любых двух постоянных значениях ${displaystyle k_ {x}}$ и ${displaystyle k_ {y},}$ это можно написать:

${displaystyle operatorname {Cov} (X, Y) = operatorname {Cov} (X-k_ {x}, Y-k_ {y}) = displaystyle {frac {sum _ {i = 1} ^ {n} (x_ { i} -k_ {x}) (y_ {i} -k_ {y}) - (сумма _ {i = 1} ^ {n} (x_ {i} -k_ {x})) (сумма _ {i = 1} ^ {n} (y_ {i} -k_ {y})) / n} {n}}.!}$

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

def shift_data_covariance(data_x, data_y):    п = len(data_x)    если п < 2:        возвращаться 0    kx = data_x[0]    ты = data_y[0]    Бывший = Эй = Exy = 0    за ix, иу в застегивать(data_x, data_y):        Бывший += ix - kx        Эй += иу - ты        Exy += (ix - kx) * (иу - ты)    возвращаться (Exy - Бывший * Эй / п) / п

Двухходовой

Двухпроходный алгоритм сначала вычисляет выборочные средние, а затем ковариацию:

${displaystyle {ar {x}} = сумма displaystyle _ {i = 1} ^ {n} x_ {i} / n}$

${displaystyle {ar {y}} = сумма displaystyle _ {i = 1} ^ {n} y_ {i} / n}$

${displaystyle operatorname {Cov} (X, Y) = displaystyle {frac {sum _ {i = 1} ^ {n} (x_ {i} - {ar {x}}) (y_ {i} - {ar {y) }})} {n}}.!}$

Двухпроходный алгоритм можно записать как:

def two_pass_covariance(data1, данные2):    п = len(data1)    среднее1 = сумма(data1) / п    среднее2 = сумма(данные2) / п    ковариация = 0    за i1, i2 в застегивать(data1, данные2):        а = i1 - среднее1        б = i2 - среднее2        ковариация += а * б / п    возвращаться ковариация

Немного более точная версия с компенсацией выполняет полный наивный алгоритм для остатков. Окончательные суммы ${displaystyle extstyle sum x_ {i}}$ и ${displaystyle extstyle sum y_ {i}}$ должен быть нулевым, но второй проход компенсирует любую небольшую ошибку.

В сети

Существует стабильный однопроходный алгоритм, похожий на онлайн-алгоритм для вычисления дисперсии, который вычисляет совместный момент ${displaystyle extstyle C_ {n} = sum _ {i = 1} ^ {n} (x_ {i} - {ar {x}} _ {n}) (y_ {i} - {ar {y}} _ { n})}$:

${displaystyle {egin {alignat} {2} {ar {x}} _ {n} & = {ar {x}} _ {n-1} &, +, & {frac {x_ {n} - {ar { x}} _ {n-1}} {n}} [5pt] {ar {y}} _ {n} & = {ar {y}} _ {n-1} &, +, & {frac { y_ {n} - {ar {y}} _ {n-1}} {n}} [5pt] C_ {n} & = C_ {n-1} &, +, & (x_ {n} - { ar {x}} _ {n}) (y_ {n} - {ar {y}} _ {n-1}) [5pt] & = C_ {n-1} &, +, & (x_ {n } - {ar {x}} _ {n-1}) (y_ {n} - {ar {y}} _ {n}) end {alignat}}}$

Кажущаяся асимметрия в этом последнем уравнении связана с тем, что ${displaystyle extstyle (x_ {n} - {ar {x}} _ {n}) = {frac {n-1} {n}} (x_ {n} - {ar {x}} _ {n-1} )}$, поэтому оба условия обновления равны ${displaystyle extstyle {frac {n-1} {n}} (x_ {n} - {ar {x}} _ {n-1}) (y_ {n} - {ar {y}} _ {n-1 })}$. Еще большей точности можно достичь, сначала вычислив средние, а затем используя устойчивый однопроходный алгоритм для остатков.

Таким образом, ковариацию можно вычислить как

${displaystyle {egin {align} operatorname {Cov} _ {N} (X, Y) = {frac {C_ {N}} {N}} & = {frac {operatorname {Cov} _ {N-1} (X , Y) cdot (N-1) + (x_ {n} - {ar {x}} _ {n}) (y_ {n} - {ar {y}} _ {n-1})} {N} } & = {frac {operatorname {Cov} _ {N-1} (X, Y) cdot (N-1) + (x_ {n} - {ar {x}} _ {n-1}) (y_ {n} - {ar {y}} _ {n})} {N}} & = {frac {operatorname {Cov} _ {N-1} (X, Y) cdot (N-1) + {frac {N-1} {N}} (x_ {n} - {ar {x}} _ {n-1}) (y_ {n} - {ar {y}} _ {n-1})} {N }} & = {frac {operatorname {Cov} _ {N-1} (X, Y) cdot (N-1) + {frac {N} {N-1}} (x_ {n} - {ar { x}} _ {n}) (y_ {n} - {ar {y}} _ {n})} {N}}. конец {выровнено}}}$

def online_covariance(data1, данные2):    среднее значение = много = C = п = 0    за Икс, у в застегивать(data1, данные2):        п += 1        dx = Икс - среднее значение        среднее значение += dx / п        много += (у - много) / п        C += dx * (у - много)    Population_covar = C / п    # Поправка Бесселя на дисперсию выборки    sample_covar = C / (п - 1)

Небольшая модификация также может быть сделана для вычисления взвешенной ковариации:

def online_weighted_covariance(data1, данные2, data3):    среднее значение = много = 0    wsum = wsum2 = 0    C = 0    за Икс, у, ш в застегивать(data1, данные2, data3):        wsum += ш        wsum2 += ш * ш        dx = Икс - среднее значение        среднее значение += (ш / wsum) * dx        много += (ш / wsum) * (у - много)        C += ш * dx * (у - много)    Population_covar = C / wsum    # Поправка Бесселя на дисперсию выборки    # Частотные веса    sample_frequency_covar = C / (wsum - 1)    # Веса надежности    sample_reliability_covar = C / (wsum - wsum2 / wsum)

Точно так же есть формула для объединения ковариаций двух наборов, которую можно использовать для распараллеливания вычислений:^[3]

${displaystyle C_ {X} = C_ {A} + C_ {B} + ({ar {x}} _ {A} - {ar {x}} _ {B}) ({ar {y}} _ {A } - {ar {y}} _ {B}) cdot {frac {n_ {A} n_ {B}} {n_ {X}}}.}$

Пакетная взвешенная версия

Также существует версия взвешенного онлайн-алгоритма с пакетным обновлением: пусть${displaystyle w_ {1}, точки w_ {N}}$ обозначим веса и напишем

${displaystyle {egin {alignat} {2} {ar {x}} _ {n + k} & = {ar {x}} _ {n} &, +, & {frac {sum _ {i = n + 1) } ^ {n + k} w_ {i} (x_ {i} - {ar {x}} _ {n})} {сумма _ {i = 1} ^ {n + k} w_ {i}}} {ar {y}} _ {n + k} & = {ar {y}} _ {n} &, +, & {frac {sum _ {i = n + 1} ^ {n + k} w_ {i } (y_ {i} - {ar {y}} _ {n})} {сумма _ {i = 1} ^ {n + k} w_ {i}}} C_ {n + k} & = C_ { n} &, +, & sum _ {i = n + 1} ^ {n + k} w_ {i} (x_ {i} - {ar {x}} _ {n + k}) (y_ {i} - {ar {y}} _ {n}) & = C_ {n} &, +, & sum _ {i = n + 1} ^ {n + k} w_ {i} (x_ {i} - {ar { x}} _ {n}) (y_ {i} - {ar {y}} _ {n + k}) end {alignat}}}$

Затем ковариацию можно вычислить как

${displaystyle operatorname {Cov} _ {N} (X, Y) = {frac {C_ {N}} {sum _ {i = 1} ^ {N} w_ {i}}}}$

Смотрите также

• Алгоритм суммирования Кахана
• Квадратные отклонения от среднего
• Ямартино метод

Рекомендации

внешняя ссылка

• Вайсштейн, Эрик В. «Пример расчета дисперсии». MathWorld.