如何计算累积几何平均值？

Question

我正在寻找使用 O(1) 空间计算累积几何平均值，其中输入是未知数字流。 Batch = 1。我有这个例子，它与我试图解决的问题足够接近：

count = 1
geo_avg = 0
while True:
  number = input("Enter number: ")

  geo_avg = ((1 + float(number)) * (1 + geo_avg)) ** (1 / count) - 1

  print("current geometric average", geo_avg, "\n")

  count += 1

我知道这是错误的。我想我可能必须删除以前的几何平均值，然后计算当前的几何平均值，但这听起来也不对。关于什么是正确使用的方程式的任何评论/提示都会非常有用。

Answer 1

您可以跟踪对数和和数量，以便将几何平均值计算为对数刻度的算术平均值，如维基百科中所述。

这是使用 clojure 的 Python 实现。

from math import log, exp


def make_geometric_averager():
    log_sum = 0
    number_count = 0

    def geometric_mean(x: float):
        nonlocal log_sum, number_count
        log_sum += log(x)
        number_count += 1
        return exp(log_sum / number_count)
    
    return geometric_mean


# Example
geometric_averager = make_geometric_averager()

print(geometric_averager(10))  # 10.0
print(geometric_averager(100))  # 31.62
print(geometric_averager(1000))  # 99.99, numerical imprecision

Answer 2

扩展评论中的想法：

import math


def running_geo_mean(x):
    sum_log = 0
    count = 0
    mean = []
    for elem in x:
        sum_log += math.log(elem)
        count += 1
        mean.append(math.exp(sum_log / count))
    return mean

使用示例：

x = [i + 1 for i in range(10)]
print(running_geo_mean(x))

[1.0, 1.414213562373095, 1.8171205928321397, 2.213363839400643, 2.6051710846973517, 2.993795165523909, 3.3800151591412964, 3.764350599503129, 4.1471662743969135, 4.528728688116766]

Answer 3

如果 g_n 是 n 个数字的几何平均值，那么包含另一个数字 x_{n+1}，新的几何平均值是

g_{n+1} = (x_{n+1} * ((g_n) ^ n)) ^ (1/(n+1))

这个公式在实践中可能会引入很多舍入误差。如果您保留正在运行的产品

(g_n) ^ n = prod{i = 1 to n} x_i = P_n

它应该返回更好的准确性。如果要包含 m 个数字，x_{n+1},...,x_{n+m}，然后计算

p = prod_{i = n+1 to n+m} x_i

先然后

P_{n+m} = (p * P_n)

g_{n+m} = P_{n+m} ^ (1/(n+m))

如果存在 P_n 溢出的危险，那么您可以按照他的建议对日志求和。我预计这会更慢，并且通常具有更高的舍入误差，但溢出的可能性要小得多。

L_n = sum_{i = 1 to n} ln(x_i) = n * ln(g_n)

g_n = exp(L_n/n)

然后再计算m个数

l = sum_{i = n+1 to n+m} ln(x_i)

L_{n+m} = l + L_n

g_{n+m} = exp(L_{n+m}/(n+m))

Answer 4

要计算大小

n

的运行几何平均值，您可以跟踪最后一个

n

数字和当前产品。当你有一个新数字时，除掉要替换的数字，然后乘以要添加的数字。

这是在 Python 中执行此操作的一种方法。

import functools
import math


def produce_geometric_average(n: int):
    lst = []
    log_prod = 0
    idx = 0

    def geometric_average(x: float):
        nonlocal lst
        nonlocal log_prod
        nonlocal idx
        log_x = math.log(x)
        if len(lst) < n:
            lst.append(log_x)
            log_prod = functools.reduce(lambda x, y: x + y, lst)
        else:
            idx_to_replace = idx % n
            old_log_x = lst[idx_to_replace]
            log_prod = log_prod - old_log_x + log_x
            lst[idx_to_replace] = log_x
        idx += 1
        return math.exp(log_prod / len(lst))
        # c.f., return math.exp(functools.reduce(lambda x, y: x + y, lst) / len(lst))

    return geometric_average


geometric_average = produce_geometric_average(3)

for x in [1.0, 1.0, 8.0, 27.0, 64.0, 125.0, 216.0]:
    print(geometric_average(x))

输出

1.0
1.0
2.0
6.000000000000001 
23.999999999999993
59.999999999999986
119.99999999999997

Answer 5

这通常是单行的，给定前一个元组 (gm_n, n) 和新值 x，返回新元组 (gm_n+1, n+1)

def update_geo_mean(gm_n, x):
    gm, n = gm_n
    return ((gm**n * x)**(1.0/(n+1)), n+1)

可以使用 exp/log 来避免上溢/下溢。

从 (1,0) 元组开始

带有日志/exp的版本

import numpy as np

def update_geo_mean(gm_n, x):
    gm, n = gm_n
    return (np.exp((n * np.log(gm) + np.log(x))/(n+1)), n+1)