python - 如何测试我的代码的数值稳定性？

Question

我有以下计算数字流的平均值和标准差的方法

import numpy as np

class RunningStatisticsVar:
    def __init__(self, ddof=0):
        self.mean = 0
        self.var = 0
        self.std = 0

        self._n = 0
        self._s = 0
        self._ddof = ddof

    def update(self, values):
        values = np.array(values, ndmin=1)
        n = len(values)

        self._n += n

        old_mean = self.mean
        delta = values - self.mean
        self.mean += (delta / self._n).sum()

        self._s += (delta * (values - self.mean)).sum()
        self.var = self._s / (self._n - self._ddof) if self._n > self._ddof else 0
        self.std = np.sqrt(self.var)

    def update_single(self, value):
        self._n += 1

        old_mean = self.mean
        self.mean += (value - old_mean) / self._n

        self._s += (value - old_mean) * (value - self.mean)
        self.var = self._s / (self._n - self._ddof) if self._n > self._ddof else 0
        self.std = np.sqrt(self.var)

    def __str__(self):
        if self.std:
            return f"(\u03BC \u00B1 \u03C3): {self.mean} \u00B1 {self.std}"
        else:
            return f"{self.mean}"

因为我想在一般应用程序中使用这个类，所以我想让它尽可能健壮。我已尽我所能通过尽可能晚地放置求和和乘法（即可能破坏数字的操作）来使代码在数值上稳定。

我有以下代码，我编写该代码是为了测试我update从 中执行批量更新的正确性update_single，我知道这是正确的：

samples = np.random.uniform(0, 100, [100000]) #
s1 = RunningStatisticsVar()
s2 = RunningStatisticsVar()

s1.update(samples)
print(s1)
for i in samples:
    s2.update_single(i)
print(s2)

为了尝试测试它的数值稳定性，我尝试使用一些非常小和一些非常大的数字，方法是将#-marked 行替换为：

samples = np.hstack((np.random.uniform(0, 1e-20, [100000]), np.random.uniform(1e20, 1e21, [100000]), np.random.uniform(0, 1e-20, [100000])))

这段代码打印出来：

(μ ± σ): 1.8363120959956487e+20 ± 3.000290391225753e+20
(μ ± σ): 1.836312095995585e+20 ± 3.000290391225752e+20

如您所见，这些数字非常相似，但并不完全相同。我的代码在数值上是否稳定，有没有比我上面的更好的方法来测试它？

标签： pythonpython-3.xnumpynumerical-stability