Блог

Я попытался смоделировать случайное число с помощью numpy.random.normal . Из этого случайного числа (среднее значение=0, std=1)

1. Я рисую несколько образцов одинакового размера (например, m=100)
2. Я вычисляю std каждого образца
3. Я беру среднее значение всех стандартных отклонений

Теоретическая статистика, а также R говорит мне, что это должно сходиться к выбранному std (который равен 1). Но почему-то, используя numpy (и scipy.stats), это не так.

Этот код генерирует рисунок, показывающий это странное поведение:

 import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, tstd

# system setup
m = 100         # number of measurments
sigma = 1       # sensor std

ez = np.arange(1,6,.05)
sample_sizes = [int(10**e) for e in ez]

# testing normal and std - they seem to work fine
sig_est = []
for n in sample_sizes:
    sample = np.random.normal(0, sigma, (n*m))
    sig_est  = [np.std(sample)]
plt.plot(ez, sig_est, marker='.', color='b', ls='', label='numpy - no means')

# numpy implementation of problem
sig_est = []
for n in sample_sizes:
    sample = np.random.normal(0, sigma, (n,m))
    sigma_est = np.std(sample, axis=1)
    sig_est  = [np.mean(sigma_est)]
plt.plot(ez, sig_est, marker='.', color='k', ls='', label='numpy')

# scipy.stats implementation
sig_est = []
for n in sample_sizes:
    sample = norm.rvs(loc=0, scale=sigma, size=(n,m))
    sigma_est = tstd(sample, axis=1)
    sig_est  = [np.mean(sigma_est)]
plt.plot(ez, sig_est, marker='.', color='r', ls='', label='scipy.stats')

plt.gca().set(xlabel = 'Number of samples [log10]')
plt.gca().legend()
plt.gca().grid(color='.9')
plt.show()
 

Есть какие-нибудь идеи?

Это интересная проблема, потому что это не проблема генератора случайных чисел, а математическая 🙂 Короткий ответ заключается в том, что все работает так, как ожидалось.

Дело в том, что в первом примере вы берете все большую и большую выборку гауссианцев i.i.d и вычисляете их стандартное отклонение с использованием np.std . Это сходится к 1, как показывает ваш график.

На втором графике вы вычисляете стандартное отклонение, всегда превышающее 100 элементов, а затем усредняете их. Таким образом, вы вычисляете не предельное значение std по многим элементам, а смещение оценки стандартного отклонения. Что, как вы выяснили, не равно нулю! Это происходит по двум причинам:

• Реализация стандартного отклонения numpy по умолчанию представляет собой квадратный корень из оценки дисперсии, которая минимизирует квадратичный риск (т. е. 1/n сумма квадратичных ошибок). Это не объективная оценка дисперсии, которая начиналась бы с 1/(n-1). Вы можете получить последнее , передав параметр ddof=1 np.std , см. Документацию здесь: https://numpy.org/doc/stable/reference/generated/numpy.std.html.
• … но даже если бы вы это сделали, вы бы не получили 0 смещений. Это потому, что вы строите график std, а не дисперсии; то есть, чтобы получить ровно 1, вы должны выровнять результаты после вычисления np.std и перед тем, как принимать среднее значение. Вы можете увидеть это, если замените свою строку

 sig_est  = [np.mean(sigma_est)]  # equivalent to sig_est.append(np.mean(sigma_est))
 

Автор:

 sig_est.append(np.mean(np.std(sample, axis=1, ddof=1)**2))
 

во втором блоке вашего кода вы действительно получите сходимость к 1.

Что касается последней реализации с использованием scipy, то, похоже, в ней используется еще одна нормализация: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.tstd.html

Они называют это «непредвзятостью», но это явно не так, с одной стороны, потому что ваши графики четко показывают это, с другой стороны, потому что точный коэффициент для получения непредвзятой оценки (для гауссов) намного сложнее, чем n/(n-1), см. Здесь: https://en.wikipedia.org/wiki/Unbiased_estimation_of_standard_deviation