Алгоритмы расчета дисперсии

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

Наивный алгоритм [ править ]

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

${\ displaystyle \ sigma ^ {2} = {\ overline {(x ^ {2})}} - {\ bar {x}} ^ {2} = \ displaystyle {\ frac {\ sum _ {i = 1} ^ {N} x_ {i} ^ {2} - (\ sum _ {i = 1} ^ {N} x_ {i}) ^ {2} / N} {N}}. \!}$

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

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

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

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

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

Вычисление сдвинутых данных [ править ]

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

${\ displaystyle \ operatorname {Var} (XK) = \ 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}}. \!}$

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

Если взять только первый образец, алгоритм можно записать на языке программирования Python как${\ displaystyle K}$

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

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

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

К  =  п  =  Ex  =  Ex2  =  0,0def  add_variable ( x ):  глобальный  K ,  n ,  Ex ,  Ex2,  если  n  ==  0 :  K  =  x  n  + =  1  Ex  + =  x  -  K  Ex2  + =  ( x  -  K )  *  ( x  -  K )def  remove_variable ( x ):  глобальный  K ,  n ,  Ex ,  Ex2  n  - =  1  Ex  - =  x  -  K  Ex2  - =  ( x  -  K )  *  ( x  -  K )def  get_mean ():  global  K ,  n ,  Ex  return  K  +  Ex  /  ndef  get_variance ():  global  n ,  Ex ,  Ex2  return  ( Ex2  -  ( Ex  *  Ex )  /  n )  /  ( n  -  1 )

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

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

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

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

${\ displaystyle {\ text {образец дисперсии}} = s ^ {2} = \ displaystyle {\ frac {\ sum _ {i = 1} ^ {n} (x_ {i} - {\ bar {x}}) ^ {2}} {n-1}}, \!}$

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

def  two_pass_variance ( данные ):  n  =  sum1  =  sum2  =  0 для  x  в  данных :  n  + =  1  sum1  + =  x среднее  =  сумма1  /  п для  x  в  данных :  sum2  + =  ( x  -  среднее )  *  ( x  -  среднее ) дисперсия  =  сумма2  /  ( n  -  1 )  дисперсия возврата 

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

Онлайн-алгоритм Велфорда [ править ]

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

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

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

${\displaystyle s_{n}^{2}={\frac {n-2}{n-1}}\,s_{n-1}^{2}+{\frac {(x_{n}-{\bar {x}}_{n-1})^{2}}{n}}=s_{n-1}^{2}+{\frac {(x_{n}-{\bar {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}-{\bar {x}}_{n-1})(x_{n}-{\bar {x}}_{n})}{n}}=\sigma _{n-1}^{2}+{\frac {(x_{n}-{\bar {x}}_{n-1})(x_{n}-{\bar {x}}_{n})-\sigma _{n-1}^{2}}{n}}.}$

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

${\displaystyle {\begin{aligned}M_{2,n}&=M_{2,n-1}+(x_{n}-{\bar {x}}_{n-1})(x_{n}-{\bar {x}}_{n})\\[4pt]s_{n}^{2}&={\frac {M_{2,n}}{n-1}}\\[4pt]\sigma _{n}^{2}&={\frac {M_{2,n}}{n}}\end{aligned}}}$

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

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

# Для нового значения newValue вычислить новый счетчик, новое среднее значение, новое значение M2. # mean накапливает среднее значение всего набора данных # M2 объединяет квадрат расстояния от среднего # count объединяет количество образцов, замеченных на данный момент def  update ( existingAggregate ,  newValue ):  ( count ,  mean ,  M2 )  =  existingAggregate  count  + =  1  delta  =  newValue  -  среднее  значение  + =  delta  /  count  delta2  =  newValue  -  среднее значение M2  + =  delta  *  delta2  return  ( количество ,  среднее ,  M2 )# Получить среднее значение, дисперсию и дисперсию выборки из агрегата def  finalize ( existingAggregate ):  ( count ,  mean ,  M2 )  =  existingAggregate  if  count  <  2 :  return  float ( "nan" )  else :  ( mean ,  variance ,  sampleVariance )  =  ( среднее ,  M2  /  число ,  M2  /  ( число  -  1 )) return  ( среднее ,  дисперсия ,  sampleVariance )

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

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

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

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

def  weighted_incremental_variance ( data_weight_pairs ):  w_sum  =  w_sum2  =  среднее значение  =  S  =  0 for  x ,  w  в  data_weight_pairs :  # В качестве альтернативы "for x, w in zip (data, weights):"  w_sum  =  w_sum  +  w  w_sum2  =  w_sum2  +  w  *  w  mean_old  =  mean  mean  =  mean_old  +  ( w  /  w_sum )  *  ( x  -  среднее_старое )  S  =  S  +  w  *  ( x  -  среднее_старое )  *  ( x -  значит ) 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 {\begin{aligned}n_{AB}&=n_{A}+n_{B}\\\delta &={\bar {x}}_{B}-{\bar {x}}_{A}\\{\bar {x}}_{AB}&={\bar {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{aligned}}}$.

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

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

def  parallel_variance ( n_a ,  avg_a ,  M2_a ,  n_b ,  avg_b ,  M2_b ):  n  =  n_a  +  n_b  delta  =  avg_b  -  avg_a  M2  =  M2_a  +  M2_b  +  delta  **  2  *  n_a  *  n_b  /  n  var_ab  =  M2  /  ( n  -  1 )  return  var_ab

Это можно обобщить, чтобы обеспечить распараллеливание с AVX , с графическими процессорами и компьютерными кластерами , а также на ковариацию. ^[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. Это серьезная проблема с наивным алгоритмом, которая возникает из-за катастрофической отмены при вычитании двух одинаковых чисел на заключительном этапе алгоритма.

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

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

${\displaystyle {\begin{aligned}M_{3,X}=M_{3,A}+M_{3,B}&{}+\delta ^{3}{\frac {n_{A}n_{B}(n_{A}-n_{B})}{n_{X}^{2}}}+3\delta {\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}\right)}{n_{X}^{3}}}\\[6pt]&{}+6\delta ^{2}{\frac {n_{A}^{2}M_{2,B}+n_{B}^{2}M_{2,A}}{n_{X}^{2}}}+4\delta {\frac {n_{A}M_{3,B}-n_{B}M_{3,A}}{n_{X}}}\end{aligned}}}$

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

${\displaystyle {\begin{aligned}&{\text{skewness}}=g_{1}={\frac {{\sqrt {n}}M_{3}}{M_{2}^{3/2}}},\\[4pt]&{\text{kurtosis}}=g_{2}={\frac {nM_{4}}{M_{2}^{2}}}-3.\end{aligned}}}$

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

${\displaystyle {\begin{aligned}\delta &=x-m\\[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 {3\delta M_{2}}{n}}\\[5pt]M_{4}'&=M_{4}+{\frac {\delta ^{4}(n-1)(n^{2}-3n+3)}{n^{3}}}+{\frac {6\delta ^{2}M_{2}}{n^{2}}}-{\frac {4\delta M_{3}}{n}}\end{aligned}}}$

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

Пример онлайн-алгоритма эксцесса, реализованного, как описано, выглядит следующим образом:

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

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

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

${\displaystyle H(x_{k})={\frac {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})\,\Delta x_{k}={\frac {1}{A}}\sum _{k=1}^{K}x_{k}^{n}h(x_{k})\,\Delta x_{k}}$

${\displaystyle \theta _{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)}{\Big )}^{n}h(x_{k})\,\Delta 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 \theta _{n}^{(h)}={\frac {1}{I}}\sum _{k=1}^{K}{\Big (}x_{k}-m_{1}^{(h)}{\Big )}^{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 {\textrm {for}}\quad n=1,2,3,4\quad {\text{ and }}\quad q=1,2,\dots ,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 {\text{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 {\text{for }}n=1,2,3,4}$

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

${\displaystyle \mu _{c}=m_{1,c}\qquad \sigma _{c}^{2}=\theta _{2,c}\qquad \alpha _{3,c}={\frac {\theta _{3,c}}{\sigma _{c}^{3}}}\qquad \alpha _{4,c}={\frac {\theta _{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  naive_covariance ( данные1 ,  данные2 ):  n  =  len ( данные1 )  sum12  =  0  sum1  =  sum ( data1 )  sum2  =  sum ( data2 ) для  i1 ,  i2  в  zip ( data1 ,  data2 ):  sum12  + =  i1  *  i2 ковариация  =  ( sum12  -  sum1  *  sum2  /  п )  /  п  возвращение  ковариационная

С оценкой среднего [ править ]

Что же касается дисперсии, ковариация двух случайных величин также инвариантное относительно сдвига, так что с учетом любых двух значений постоянных и его можно записать:${\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})-(\sum _{i=1}^{n}(x_{i}-k_{x}))(\sum _{i=1}^{n}(y_{i}-k_{y}))/n}{n}}.\!}$

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

def  shift_data_covariance ( data_x ,  data_y ):  n  =  len ( data_x )  if  n  <  2 :  return  0  kx  =  data_x [ 0 ]  ky  =  data_y [ 0 ]  Ex  =  Ey  =  Exy  =  0  for  ix ,  iy  in  zip ( data_x ,  data_y ):  Ex  + =  ix  -  kx Ey  + =  iy  -  ky  Exy  + =  ( ix  -  kx )  *  ( iy  -  ky )  return  ( Exy  -  Ex  *  Ey  /  n )  /  n

Двухходовой [ править ]

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

${\displaystyle {\bar {x}}=\displaystyle \sum _{i=1}^{n}x_{i}/n}$

${\displaystyle {\bar {y}}=\displaystyle \sum _{i=1}^{n}y_{i}/n}$

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

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

def  two_pass_covariance ( данные1 ,  данные2 ):  n  =  len ( данные1 ) среднее1  =  сумма ( данные1 )  /  п  среднее2  =  сумма ( данные2 )  /  п ковариация  =  0 для  i1 ,  i2  в  zip ( data1 ,  data2 ):  a  =  i1  -  mean1  b  =  i2  -  mean2  covariance  + =  a  *  b  /  n  return  covariance

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

Онлайн [ править ]

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

${\displaystyle {\begin{alignedat}{2}{\bar {x}}_{n}&={\bar {x}}_{n-1}&\,+\,&{\frac {x_{n}-{\bar {x}}_{n-1}}{n}}\\[5pt]{\bar {y}}_{n}&={\bar {y}}_{n-1}&\,+\,&{\frac {y_{n}-{\bar {y}}_{n-1}}{n}}\\[5pt]C_{n}&=C_{n-1}&\,+\,&(x_{n}-{\bar {x}}_{n})(y_{n}-{\bar {y}}_{n-1})\\[5pt]&=C_{n-1}&\,+\,&(x_{n}-{\bar {x}}_{n-1})(y_{n}-{\bar {y}}_{n})\end{alignedat}}}$

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

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

${\displaystyle {\begin{aligned}\operatorname {Cov} _{N}(X,Y)={\frac {C_{N}}{N}}&={\frac {\operatorname {Cov} _{N-1}(X,Y)\cdot (N-1)+(x_{n}-{\bar {x}}_{n})(y_{n}-{\bar {y}}_{n-1})}{N}}\\&={\frac {\operatorname {Cov} _{N-1}(X,Y)\cdot (N-1)+(x_{n}-{\bar {x}}_{n-1})(y_{n}-{\bar {y}}_{n})}{N}}\\&={\frac {\operatorname {Cov} _{N-1}(X,Y)\cdot (N-1)+{\frac {N-1}{N}}(x_{n}-{\bar {x}}_{n-1})(y_{n}-{\bar {y}}_{n-1})}{N}}\\&={\frac {\operatorname {Cov} _{N-1}(X,Y)\cdot (N-1)+{\frac {N}{N-1}}(x_{n}-{\bar {x}}_{n})(y_{n}-{\bar {y}}_{n})}{N}}.\end{aligned}}}$

Защиту  online_covariance ( data1 ,  data2 ):  meanx  =  Meany  =  С  =  п  =  0  для  х ,  у  в  молнии ( DATA1 ,  data2 ):  п  + =  1  дх  =  х  -  meanx  meanx  + =  дх  /  п  Meany  + =  ( у  -  meany )  /  n  C  + =  dx  *  ( y  - meany ) population_covar  =  С  /  п  коррекции # Бесселя для образца дисперсии  sample_covar  =  С  /  ( п  -  1 )

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

Защиту  online_weighted_covariance ( data1 ,  data2 ,  данные3 ):  meanx  =  Meany  =  0  wsum  =  wsum2  =  0  С  =  0  для  х ,  у ,  ш  в  молнии ( data1 ,  data2 ,  Data3 ):  wsum  + =  ш  wsum2  + =  ш  *  ш  дх  =  x  -  среднее  x среднее x  + =  ( w /  wsum )  *  dx  meany  + =  ( w  /  wsum )  *  ( y  -  средний )  C  + =  w  *  dx  *  ( y  -  средний ) population_covar  =  С  /  wsum  # Бесселя коррекции для образца дисперсии  # частот весов  sample_frequency_covar  =  С  /  ( wsum  -  1 )  # Надежность веса  sample_reliability_covar  =  С  /  ( wsum  -  wsum2  /  wsum )

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

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

Пакетная взвешенная версия [ править ]

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

${\displaystyle {\begin{alignedat}{2}{\bar {x}}_{n+k}&={\bar {x}}_{n}&\,+\,&{\frac {\sum _{i=n+1}^{n+k}w_{i}(x_{i}-{\bar {x}}_{n})}{\sum _{i=1}^{n+k}w_{i}}}\\{\bar {y}}_{n+k}&={\bar {y}}_{n}&\,+\,&{\frac {\sum _{i=n+1}^{n+k}w_{i}(y_{i}-{\bar {y}}_{n})}{\sum _{i=1}^{n+k}w_{i}}}\\C_{n+k}&=C_{n}&\,+\,&\sum _{i=n+1}^{n+k}w_{i}(x_{i}-{\bar {x}}_{n+k})(y_{i}-{\bar {y}}_{n})\\&=C_{n}&\,+\,&\sum _{i=n+1}^{n+k}w_{i}(x_{i}-{\bar {x}}_{n})(y_{i}-{\bar {y}}_{n+k})\\\end{alignedat}}}$

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

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

См. Также [ править ]

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

Ссылки [ править ]

Внешние ссылки [ править ]

• Вайсштейн, Эрик В. "Расчет дисперсии выборки" . MathWorld .