TL; DR
Логнормальное распределение с положительным перекосом и распределением с тяжелыми хвостами.При выполнении арифметических операций с плавающей запятой (таких как сумма, среднее или стандартное отклонение) для выборки, взятой из сильно искаженного распределения, вектор выборки содержит значения с расхождением в несколько порядков (много десятилетий).Это делает вычисления неточными.
Проблема возникает из этих двух строк:
mean_calc += [np.average(ln)]
sigma_calc += [np.std(ln)]
Поскольку ln содержит как очень низкие, так и очень высокие значения с порядком величины, намного превышающим точность с плавающей запятой.
Проблема может быть легко обнаружена, чтобы предупредить пользователя, что ее вычисления неверны, используя следующий предикат:
(max(ln) + min(ln)) <= max(ln)
Что явно неверно в Strictly Positive Real, но должно учитываться при использовании FiniteПрецизионная арифметика.
Изменение вашего MCVE
Если мы слегка изменим ваш MCVE на:
from scipy import stats
for s in ss:
mu = -0.5*s*s
ln = stats.lognorm(s, scale=np.exp(mu)).rvs(N*N)
f = stats.lognorm.fit(ln, floc=0)
mean_calc += [f[2]*np.exp(0.5*s*s)]
sigma_calc += [np.sqrt((np.exp(f[0]**2)-1)*(np.exp(2*mu + s*s)))]
mean_analytic += [np.exp(mu+0.5*s*s)]
sigma_analytic += [np.sqrt((np.exp(s**2)-1)*(np.exp(2*mu + s*s)))]
Это дает достаточно правильное среднее значение и стандартное отклонениеоценка даже для высокого значения сигмы.
Ключ в том, что fit использует алгоритм MLE для оценки параметров.Это полностью отличается от вашего первоначального подхода, который непосредственно выполняет среднее значение выборки.
Метод fit возвращает кортеж с (sigma, loc=0, scale=exp(mu)), которые являются параметрами объекта scipy.stats.lognorm какуказано в документации.
Я думаю, вам следует выяснить, как вы оцениваете среднее значение и стандартное отклонение.Расхождение, вероятно, происходит из этой части вашего алгоритма.
Может быть несколько причин, по которым оно расходится, по крайней мере, рассмотрим:
- Смещенная оценка: Является ли ваша оценка правильной инесмещенный? Среднее - это несмещенная оценка (см. следующий раздел), но, возможно, неэффективная;
- Выбранные значения от псевдослучайного генератора могут быть не такими интенсивными, как их следует сравнивать с теоретическим распределением: возможноMLE менее чувствителен, чем ваша оценка Новый сильфон MCVE не поддерживает эту гипотезу, но арифметическая ошибка с плавающей точкой может объяснить, почему ваши оценки занижены;
- арифметическая ошибка с плавающей точкой Новая сильфонная подсветка MCVEчто это является частью вашей проблемы.
Научная цитата
Кажется, наивная оценка среднего значения (просто принимая среднее значение), даже если она беспристрастна, неэффективна для правильной оценки среднего значения для большой сигмы(см. Ци Тан , Сравнение различных методов для оценкиОбычные средства , с.11):
Наивный оценщик легко рассчитать и он беспристрастен.Тем не менее, эта оценка может быть неэффективной, когда дисперсия велика и размер выборки невелик.
В диссертации рассматривается несколько методов оценки среднего логнормального распределения и используется MLE в качестве эталона для сравнения.Это объясняет, почему ваш метод имеет дрейф, так как сигма возрастает, а MLE лучше придерживается, увы, он не эффективен по времени для больших N. Очень интересная статья.
Статистические соображения
Вспоминая, чем:
- Логнормальное - тяжелое и длиннохвостое распределение положительно перекошено.Одним из следствий этого является то, что с ростом сигмы параметра формы растет асимметрия и асимметрия, а также возрастает сила выбросов.
- Эффект размера выборки: по мере роста числа выборок, взятых из распределения, ожидается, чтовыбросы увеличиваются (увеличивается и экстент).
Создание нового MCVE
Позволяет построить новый MCVE, чтобы сделать его более понятным.Нижеприведенный код рисует образцы разных размеров (N в диапазоне от 100 до 10000) из логнормального распределения, где параметр формы изменяется (sigma в диапазоне от 0.1 и 10), а параметр масштаба установлен наунитарный.
import warnings
import numpy as np
from scipy import stats
# Make computation reproducible among batches:
np.random.seed(123456789)
# Parameters ranges:
sigmas = np.arange(0.1, 10.1, 0.1)
sizes = np.logspace(2, 5, 21, base=10).astype(int)
# Placeholders:
rv = np.empty((sigmas.size,), dtype=object)
xmean = np.full((3, sigmas.size, sizes.size), np.nan)
xstd = np.full((3, sigmas.size, sizes.size), np.nan)
xextent = np.full((2, sigmas.size, sizes.size), np.nan)
eps = np.finfo(np.float64).eps
# Iterate Shape Parameter:
for (i, s) in enumerate(sigmas):
# Create Random Variable:
rv[i] = stats.lognorm(s, loc=0, scale=1)
# Iterate Sample Size:
for (j, N) in enumerate(sizes):
# Draw Samples:
xs = rv[i].rvs(N)
# Sample Extent:
xextent[:,i,j] = [np.min(xs), np.max(xs)]
# Check (max(x) + min(x)) <= max(x)
if (xextent[0,i,j] + xextent[1,i,j]) - xextent[1,i,j] < eps:
warnings.warn("Potential Float Arithmetic Errors: logN(mu=%.2f, sigma=%2f).sample(%d)" % (0, s, N))
# Generate different Estimators:
# Fit Parameters using MLE:
fit = stats.lognorm.fit(xs, floc=0)
xmean[0,i,j] = fit[2]
xstd[0,i,j] = fit[0]
# Naive (Bad Estimators because of Float Arithmetic Error):
xmean[1,i,j] = np.mean(xs)*np.exp(-0.5*s**2)
xstd[1,i,j] = np.sqrt(np.log(np.std(xs)**2*np.exp(-s**2)+1))
# Log-transform:
xmean[2,i,j] = np.exp(np.mean(np.log(xs)))
xstd[2,i,j] = np.std(np.log(xs))
Наблюдение: новый MCVE начинает выдавать предупреждения, когда sigma > 4.
MLE в качестве ссылки
Оценка параметров формы и масштаба с использованием MLE дает хорошие результаты:
На двух приведенных выше рисунках показано, что:
- Ошибка при оценке растет вместе с параметром формы;
- Ошибка при оценке уменьшается с увеличением размера выборки ( CTL );
Обратите внимание, что MLE также хорошо подходит для параметра формы:
Арифметика с плавающей запятой
Стоит представить размер вытянутых образцовв зависимости от параметра формы и размера выборки:
Или десятичная величина между наименьшим и наибольшим числом в образце:
На моей установке:
np.finfo(np.float64).precision # 15
np.finfo(np.float64).eps # 2.220446049250313e-16
Это означает, что у нас есть максимум 15 значащих цифр для работы, если величина между двумя числами превышает, то наибольшее число поглощаетменьшие.
Базовый пример: Каков результат 1 + 1e6, если мы можем только сохранить fнаши значимые цифры?Точный результат - 1,000,001.0, но он должен быть округлен до 1.000e6.Это подразумевает: результат операции равен наибольшему числу из-за точности округления.Она присуща Арифметика конечной точности .
Два предыдущих рисунка в сочетании со статистическим учетом подтверждают ваше наблюдение, что увеличение N не улучшает оценку для большого значения sigma вваш MCVE.
Цифры выше и ниже показывают, что когда sigma > 3 у нас недостаточно значащих цифр (менее 5) для выполнения правильных вычислений.
Более того, мы можем сказать, что оценщик будет недооценен , так как самые большие числа будут поглощать наименьшее, а недооцененная сумма будет затем разделена на N, делая оценку смещенной по умолчанию.
Когда параметр формы становится достаточно большим, вычислениясильно смещены из-за арифметических ошибок с плавающей точкой .
Это означает использование величин, таких как:
np.mean(xs)
np.std(xs)
Когда расчетная оценка будет иметь огромную арифметическую ошибку с плавающей точкой, потому чтоважного расхождения между значениями, хранящимися в xs.Приведенные ниже цифры воспроизводят вашу проблему:
Как указано, оценки по умолчанию (не в избытке), поскольку высокиезначения (несколько выбросов) в выбранном векторе поглощают небольшие значения (большинство выборочных значений).
Логарифмическое преобразование
Если мы применяем логарифмическое преобразование , мы можем резко уменьшитьэто явление:
xmean[2,i,j] = np.exp(np.mean(np.log(xs)))
xstd[2,i,j] = np.std(np.log(xs))
И тогда наивная оценка среднего значения является правильной и будет намного меньше подвержена арифметической ошибке с плавающей запятой , поскольку все значения выборки будут лежать в течение нескольких десятилетийвместо того, чтобы иметь относительную величину, превышающую арифметическую точность с плавающей точкой .
На самом деле, лог-преобразование возвращает тот же результат для среднего и стандартного вычисления, что и MLE для каждого N и sigma:
np.allclose(xmean[0,:,:], xmean[2,:,:]) # True
np.allclose(xstd[0,:,:], xstd[2,:,:]) # True
Справочник
Если вам нужны полные и подробные объяснения такого рода проблем при выполнении научных вычислений, я рекомендую вам: NJ Higham , Точность и стабильность численных алгоритмов, Сиам, второе издание, 2002 .
Bonus
Вот пример кода генерации фигуры:
import matplotlib.pyplot as plt
fig, axe = plt.subplots()
idx = slice(None, None, 5)
axe.loglog(sigmas, xmean[0,:,idx])
axe.axhline(1, linestyle=':', color='k')
axe.set_title(r"MLE: $x \sim \log\mathcal{N}(\mu=0,\sigma)$")
axe.set_xlabel(r"Standard Deviation, $\sigma$")
axe.set_ylabel(r"Mean Estimation, $\hat{\mu}$")
axe.set_ylim([0.1,10])
lgd = axe.legend([r"$N = %d$" % s for s in sizes[idx]] + ['Exact'], bbox_to_anchor=(1,1), loc='upper left')
axe.grid(which='both')
fig.savefig('Lognorm_MLE_Emean_Sigma.png', dpi=120, bbox_extra_artists=(lgd,), bbox_inches='tight')