Как мне проверить мой код на числовую стабильность?

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)

Чтобы попытаться проверить его числовую стабильность, я попытался использовать некоторыеочень маленькие и некоторые очень большие цифры, заменив строку, помеченную #, на:

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

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