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


Из Википедии, свободной энциклопедии

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

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

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

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

Поэтому, наивный алгоритм расчета оцененной дисперсии даются ниже:

  • Пусть п ← 0, Sum ← 0, СУММКВ ← 0
  • Для каждой точки привязки х :
    • пп + 1
    • Сумма ← Sum + х
    • СУММЫ ← СУММЫ + х х х
  • Вар = (СУММЫ - (сумма × Сумма) / п) / (п - 1)

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

Поскольку СУММЫ и (Сумма × Сумма) / п могут быть очень похожие номера, отмена может привести к точности результата , чтобы быть значительно меньше , чем присущая точность арифметики с плавающей точкой , используемой для выполнения вычислений. Таким образом , этот алгоритм не следует применять на практике, и несколько альтернативных, численно стабильными, были предложены алгоритмы. Это особенно плохо , если стандартное отклонение мало по сравнению со средним значением. Тем не менее, алгоритм может быть улучшен путем применения способа в соответствии с предполагаемым средним .

Вычислительный сдвинуты данные

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

с какой - либо константой, что приводит к новой формуле

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

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

def shifted_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 - K) * (x - K)
   variance = (Ex2 - (Ex * Ex)/n)/(n - 1)
   # use n instead of (n-1) if want to compute the exact variance of the given data
   # use (n-1) if data are samples of a larger population
   return variance

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

K = n = Ex = Ex2 = 0.0

def add_variable(x):
    if (n == 0):
      K = x
    n = n + 1
    Ex += x - K
    Ex2 += (x - K) * (x - K)

def remove_variable(x):
    n = n - 1
    Ex -= (x - K)
    Ex2 -= (x - K) * (x - K)

def get_meanvalue():
    return K + Ex / n

def get_variance():
    return (Ex2 - (Ex*Ex)/n) / (n-1)

Алгоритм Двухпроходного

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

,

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

,

где s представляет собой стандартное отклонение. Это дается следующий псевдокод:

def two_pass_variance(data):
    n = sum1 = sum2 = 0

    for x in data:
        n += 1
        sum1 += x

    mean = sum1 / n

    for x in data:
        sum2 += (x - mean)*(x - mean)

    variance = sum2 / (n - 1)
    return variance

Этот алгоритм численно устойчив , если п мала. Тем не менее, результаты обоих этих простых алгоритмов ( «Наивные» и «Двухпроходное») может зависеть премного от упорядочения данных и может дать плохие результаты для очень больших наборов данных из - за повторных ошибок округления в накоплении суммы. Такие методы, как компенсированное суммирование могут быть использованы для борьбы с этой ошибкой в некоторой степень.

Интернет алгоритм Уэлфорд в

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

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

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

Этот алгоритм был найден Welford, и она была тщательно проанализирована. Кроме того , общие для обозначения и .

Реализация примера Python для алгоритма Уэлфорда будет дан ниже.

# for a new value newValue, compute the new count, new mean, the new M2.
# mean accumulates the mean of the entire dataset
# M2 aggregates the squared distance from the mean
# count aggregates the number of samples seen so far
def update(existingAggregate, newValue):
    (count, mean, M2) = existingAggregate
    count += 1 
    delta = newValue - mean
    mean += delta / count
    delta2 = newValue - mean
    M2 += delta * delta2

    return (count, mean, M2)

# retrieve the mean, variance and sample variance from an aggregate
def finalize(existingAggregate):
    (count, mean, M2) = existingAggregate
    (mean, variance, sampleVariance) = (mean, M2/count, M2/(count - 1)) 
    if count < 2:
        return float('nan')
    else:
        return (mean, variance, sampleVariance)

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

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

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

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

def weighted_incremental_variance(dataWeightPairs):
    wSum = wSum2 = mean = S = 0

    for x, w in dataWeightPairs:  # Alternatively "for x, w in zip(data, weights):"
        wSum = wSum + w
        wSum2 = wSum2 + w*w
        meanOld = mean
        mean = meanOld + (w / wSum) * (x - meanOld)
        S = S + w * (x - meanOld) * (x - mean)

    population_variance = S / wSum
    # Bessel's correction for weighted samples
    # Frequency weights
    sample_frequency_variance = S / (wSum - 1)
    # Reliability weights
    sample_reliability_variance = S / (wSum - wSum2/wSum)

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

Чен и др. обратите внимание , что выше алгоритм «Он-лайн» является частным случаем алгоритма , который работает для любого разбиения образца на множества , :

,

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

Метод Чана для оценки среднего значения численно неустойчиво , когда и как большие, так как численная ошибка не масштабируется вниз таким образом , что он находится в случае. В таких случаях предпочтение .

def parallel_variance(avg_a, count_a, var_a, avg_b, count_b, var_b):
    delta = avg_b - avg_a
    m_a = var_a * (count_a - 1)
    m_b = var_b * (count_b - 1)
    M2 = m_a + m_b + delta ** 2 * count_a * count_b / (count_a + count_b)
    return M2 / (count_a + count_b - 1)

пример

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

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

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

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

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

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

перекос:
эксцесс:

Для дополнительного случая (то есть ), это упрощает для:

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

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

def online_kurtosis(data):
    n = mean = M2 = M3 = M4 = 0

    for x in data:
        n1 = n
        n = n + 1
        delta = x - mean
        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*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

    kurtosis = (n*M4) / (M2*M2) - 3
    return kurtosis

Pébaÿ дополнительно расширяет эти результаты на произвольные порядка центральных моментов , для инкрементного и попарных случаев, а впоследствии Pébaÿ и соавт. для взвешенных и сложных моментов. Можно также найти там аналогичные формулы для ковариации .

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

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

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

Второй подход от Choi и Sweetman является аналитической методики, чтобы объединить статистические моменты из отдельных сегментов тайм-истории таким образом, что полученные в результате общие моменты являются те, полного времени истории. Эта методика может быть использована для параллельного вычисления статистических моментов с последующей комбинацией тех моментов, или для комбинации статистических моментов, вычисленных в последовательные времена.

Если известны наборы статистических моментов: для , то каждый может быть выражено в терминах эквивалентных сырых моментов:

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

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

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

Известные отношения между необработанными моментами ( ) и центральными моментами ( ), затем используются для вычисления центральных моментов каскадного времени истории. Наконец, статистические моменты каскадной истории вычисляется из центральных моментов:

ковариации

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

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

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

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

def naive_covariance(data1, data2):
    n = len(data1)
    sum12 = 0
    sum1 = sum(data1)
    sum2 = sum(data2)

    for i1, i2 in zip(data1, data2):
        sum12 += i1*i2

    covariance = (sum12 - sum1*sum2 / n) / n
    return covariance

При оценке среднего

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

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

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

Двухпроходное

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

Алгоритм два прохода может быть записан в виде:

def two_pass_covariance(data1, data2):
    n = len(data1)

    mean1 = sum(data1) / n
    mean2 = sum(data2) / n

    covariance = 0

    for i1, i2 in zip(data1, data2):
        a = i1 - mean1
        b = i2 - mean2
        covariance += a*b / n
    return covariance

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

онлайн

Устойчивый алгоритм один проход существует, похож на онлайн - алгоритм для вычисления дисперсии, которая вычисляет момент сотрудничества :

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

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

def online_covariance(data1, data2):
    meanx = meany = C = n = 0
    for x, y in zip(data1, data2):
        n += 1
        dx = x - meanx
        meanx += dx / n
        meany += (y - meany) / n
        C += dx * (y - meany)

    population_covar = C / n
    # Bessel's correction for sample variance
    sample_covar = C / (n - 1)

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

def online_weighted_covariance(data1, data2, data3):
    meanx = meany = 0
    wsum = wsum2 = 0
    C = 0
    for x, y, w in zip(data1, data2, data3):
        wsum += w
        wsum2 += w*w
        dx = x - meanx
        meanx += (w / wsum) * dx
        meany += (w / wsum) * (y - meany)
        C += w * dx * (y - meany)

    population_covar = C / wsum
    # Bessel's correction for sample variance
    # Frequency weights
    sample_frequency_covar = C / (wsum - 1)
    # Reliability weights
    sample_reliability_covar = C / (wsum - wsum2 / wsum)

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

Взвешенная порционный версия

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

Ковариация затем может быть вычислена как

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

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

  1. ^ Б Бо Einarsson (1 августа 2005). Точность и надежность в научных вычислениях . SIAM. п. 47. ISBN  978-0-89871-584-2 . Проверено 17 Февралю 2013 .
  2. ^ Б T.F.Chan, GH Голуб и RJ Левек (1983). « » Алгоритмы для вычисления дисперсии выборки: анализ и рекомендации», Американский Statistician, 37" (PDF) : 242-247.
  3. ^ Б Шуберт, Эрих; Герц, Майкл (2018-07-09). «Численно стабильно параллельное вычисление (со) дисперсии» . ACM: 10. DOI : 10,1145 / 3221269,3223036 . ISBN  9781450365055 .
  4. ^ Higham, Николай (2002). Точность и стабильность численных алгоритмов (2 изд) (проблема 1.10) . SIAM.
  5. ^ BP Welford (1962). «Обратите внимание на методе расчета скорректированных сумм квадратов и изделий» . Technometrics 4 (3): 419-420.
  6. ^ Дональда Кнута (1998). Искусство программирования , том 2: Получисленные алгоритмы ., Третьего обявили, с. 232. Boston: Addison-Wesley.
  7. ^ Чан, Тони F .; Голуб, Джин Н .; Левек, Randall J. (1983). Алгоритмы для вычисления выборочной дисперсии: анализ и рекомендация. Американский Статистик 37, 242-247. https://www.jstor.org/stable/2683386
  8. ^ Лин, Роберт Ф. (1974). Сравнение нескольких алгоритмов для вычисления выборочных средних и дисперсий. Журнал Американской статистической ассоциации, Vol. 69, № 348, 859-866. DOI : 10,2307 / 2286154
  9. ^ http://www.johndcook.com/standard_deviation.html
  10. ^ DHD West (1979). Коммуникации по АКМ , 22, 9, 532-535: Обновление среднего значения и дисперсии оценок: усовершенствовал метод
  11. ^ Чен, Тони Ф. ; Голуб, Джин H. ; Левек, Randall J. (1979), «Обновление формул и Попарный алгоритм вычисление Sample дисперсий.» (PDF) , технический отчет STAN-CS-79-773 , факультет компьютерных наук Стэнфордского университета,
  12. ^ Terriberry, Тимоти Б. (2007), Вычислительный высшего порядка Moments Интернет
  13. ^ Pébaÿ, Филипп (2008), «Формула для Robust, однопроходного параллельного вычисления ковариации и произвольных порядки статистических моментов» (PDF) , технический отчет SAND2008-6212 , Sandia National Laboratories
  14. ^ Pébaÿ, Филипп; Terriberry, Тимоти; Колла, Hemanth; Bennett, Жанин (2016), «Числено Стабильные масштабируемые формулы для параллельного и Интернет Вычисления высших порядков Многофакторных центральных моментов с произвольными весами» , Вычислительные статистики , Springer
  15. ^ Б Choi, Myoungkeun; Sweetman, Берт (2010), «Эффективное Расчет статистических моментов для мониторинга структурного здоровья» , журнал мониторинга Структурного здравоохранения , 9 (1)

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