Какова формула алгоритма Уэлфорда для дисперсии / стандартизации с пакетными обновлениями? - PullRequest
Новая
       29

Какова формула алгоритма Уэлфорда для дисперсии / стандартизации с пакетными обновлениями?

0 голосов
/ 01 июня 2019

Я хочу расширить онлайновый алгоритм Уэлфорда, чтобы его можно было обновлять несколькими номерами (в пакетном режиме) вместо одного по одному: https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance

Я пытался обновить алгоритм со страницы викикак это: 

# my attempt.
def update1(existingAggregate, newValues):
    (count, mean, M2) = existingAggregate
    count += len(newValues) 
    delta = np.sum(np.subtract(newValues, [mean] * len(newValues)))
    mean += delta / count
    delta2 = np.sum(np.subtract(newValues, [mean] * len(newValues)))
    M2 += delta * delta2

    return (count, mean, M2)

# The original two functions from wikipedia.
def update(existingAggregate, newValue):
    (count, mean, M2) = existingAggregate
    count += 1 
    delta = newValue - mean
    mean += delta / count
    delta2 = newValue - mean
    M2 += delta * delta2

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)

Однако я не должен понимать это правильно, потому что результат неправильный: 

# example x that might have led to an a = (2, 2.0, 2.0).
x = [1.0, 3.0]
mean = np.mean(x)
count = len(x)
m2 = np.sum(np.subtract(x, [mean] * count)**2)

a = (count, mean, m2)
print(a)
# new batch of values.
b = [5, 3]

Обратите внимание, что a = (2, 2.0, 2.0) означает, что мыбыло 2 наблюдения, и их среднее значение было 2,0. 

# update one at a time.
temp = update(a, newValues[0])
result_single = update(temp, newValues[1])
print(finalize(result_single))

# update with my faulty batch function.
result_batch = update1(a, newValues)
print(finalize(result_batch))

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

(3.0, 2.0, 2.6666666666666665)
(3.0, 2.5, 3.3333333333333335)

Чего мне не хватает в отношении правильногоОбновления дисперсии?Нужно ли мне как-то обновлять функцию финализации?

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

Ответы [ 2 ]

1 голос
/ 01 июня 2019

Я не слишком знаком с Python, поэтому я скорее буду придерживаться математической записи.

Чтобы обновить среднее значение, вы должны сделать: 

s = sum of new values
c = number of new values
newMean = oldMean + sum_i (newValue[i] - oldMean) / newCount

Для M2, вы должны добавить еще одно суммирование: 

newM2 = oldM2 + sum_i ((newValue[i] - newMean) * (newValue[i] - oldMean))

Я не уверен, что вы действительно сохраните что-нибудь с этим массовым обновлением, поскольку у вас все еще есть цикл внутри.

0 голосов
/ 01 июня 2019

Благодаря разъяснениям Нико Я понял это! Проблема заключалась в том, что я суммировал для дельт, а затем умножал, чтобы получить М2, но вместо этого мне пришлось суммировать по произведению дельт. Вот правильная функция пакетной обработки, которая может принимать как отдельные числа, так и пакеты: 

# https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
def update(existingAggregate, newValues):
    if isinstance(newValues, (int, float, complex)):
        # Handle single digits.
        newValues = [newValues]

    (count, mean, M2) = existingAggregate
    count += len(newValues) 
    # newvalues - oldMean
    delta = np.subtract(newValues, [mean] * len(newValues))
    mean += np.sum(delta / count)
    # newvalues - newMeant
    delta2 = np.subtract(newValues, [mean] * len(newValues))
    M2 += np.sum(delta * delta2)

    return (count, mean, M2)

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)

Пример использования: 

x = [1.0, 3.0]
mean = np.mean(x)
count = len(x)
m2 = np.sum(np.subtract(x, [mean] * count)**2)

a = (count, mean, m2)
print(a)
# new batch of values.
b = [5, 3]

result_batch = update(a, b)
result_batch1 = update(a, b[0])

print(finalize(result_batch))
print(finalize(result_batch1))

И это действительно быстрее: 

import timeit
x = random.sample(range(1, 10000), 1000)
# ...
b = random.sample(range(1, 10000), 1000)

start_time = timeit.default_timer()
result_batch = update(a, b)
print(f'{timeit.default_timer() - start_time:.4f}')
print(*(f'{x:.2f}' for x in finalize(result_batch)))

start_time = timeit.default_timer()
for i in b:
    a  = update1(a, i)
print(f'{timeit.default_timer() - start_time:.4f}')
print(*(f'{x:.2f}' for x in finalize(result_batch)))

Результат: 

0.0010
5008.36 8423224.68 8427438.40
0.0031
5008.36 8423224.68 8427438.40
...