Правильный способ получения доверительного интервала со скупым
У меня есть одномерный массив данных:
a = np.array([1,2,3,4,4,4,5,5,5,5,4,4,4,6,7,8])
для которого я хочу получить 68% доверительный интервал (т.е. 1 сигма).
Первый комментарий в этом ответе утверждает, что это может быть достигнуто с помощью scipy.stats.norm.interval из функции scipy.stats.norm через:
from scipy import stats
import numpy as np
mean, sigma = np.mean(a), np.std(a)
conf_int = stats.norm.interval(0.68, loc=mean,
scale=sigma)
Но комментарий в этом посте гласит, что фактический правильный способ получения доверительного интервала:
conf_int = stats.norm.interval(0.68, loc=mean,
scale=sigma / np.sqrt(len(a)))
то есть сигма делится на квадратный корень размера выборки: np.sqrt(len(a)),
Вопрос в том, какая версия является правильной?
3 ответа
68-процентный доверительный интервал для одиночного розыгрыша из нормального распределения со средним значением отклонения mu и std составляет
stats.norm.interval(0.68, loc=mu, scale=sigma)
68% доверительный интервал для среднего значения N взят из нормального распределения со средним значением mu и сигма стандартного отклонения
stats.norm.interval(0.68, loc=mu, scale=sigma/sqrt(N))
Интуитивно понятно, что эти формулы имеют смысл, поскольку, если вы поднимите банку желейных бобов и попросите большое количество людей угадать количество желейных бобов, у каждого человека может быть слишком много - то же самое стандартное отклонение sigma - но среднее из предположений будет делать замечательную работу по оценке фактического числа, и это отражается стандартным отклонением среднего значения, уменьшающегося с коэффициентом 1/sqrt(N),
Если один розыгрыш имеет дисперсию sigma**2 тогда по формуле Бинаеме сумма N некоррелированные ничьи имеют дисперсию N*sigma**2,
Среднее значение равно сумме, деленной на N. Когда вы умножаете случайную величину (например, сумму) на константу, дисперсия умножается на квадрат в квадрате. То есть
Var(cX) = c**2 * Var(X)
Таким образом, дисперсия среднего равна
(variance of the sum)/N**2 = N * sigma**2 / N**2 = sigma**2 / N
и поэтому стандартное отклонение от среднего значения (которое является квадратным корнем дисперсии) равно
sigma/sqrt(N).
Это происхождение sqrt(N) в знаменателе.
Вот пример кода, основанного на коде Тома, который демонстрирует утверждения, сделанные выше:
import numpy as np
from scipy import stats
N = 10000
a = np.random.normal(0, 1, N)
mean, sigma = a.mean(), a.std(ddof=1)
conf_int_a = stats.norm.interval(0.68, loc=mean, scale=sigma)
print('{:0.2%} of the single draws are in conf_int_a'
.format(((a >= conf_int_a[0]) & (a < conf_int_a[1])).sum() / float(N)))
M = 1000
b = np.random.normal(0, 1, (N, M)).mean(axis=1)
conf_int_b = stats.norm.interval(0.68, loc=0, scale=1 / np.sqrt(M))
print('{:0.2%} of the means are in conf_int_b'
.format(((b >= conf_int_b[0]) & (b < conf_int_b[1])).sum() / float(N)))
печать
68.03% of the single draws are in conf_int_a
67.78% of the means are in conf_int_b
Остерегайтесь, если вы определите conf_int_b с оценками для mean а также sigma на основании образца a среднее значение может не попасть в conf_int_b с желаемой частотой.
Если вы берете выборку из распределения и вычисляете среднее значение выборки и стандартное отклонение,
mean, sigma = a.mean(), a.std()
обратите внимание, что нет никакой гарантии, что они будут равны среднему значению и стандартному отклонению, и что мы предполагаем, что население распределено нормально - это не автоматические данные!
Если вы берете выборку и хотите оценить среднее значение по населению и стандартное отклонение, вы должны использовать
mean, sigma = a.mean(), a.std(ddof=1)
поскольку это значение для сигмы является объективной оценкой стандартного отклонения населения.
Я только что проверил, как R и GraphPad вычисляют доверительные интервалы, и они увеличивают интервал в случае небольшого размера выборки (n). Например, более чем в 6 раз для n=2 по сравнению с большим n. Этот код (на основе ответа Шасана) соответствует их доверительным интервалам:
import numpy as np, scipy.stats as st
# returns confidence interval of mean
def confIntMean(a, conf=0.95):
mean, sem, m = np.mean(a), st.sem(a), st.t.ppf((1+conf)/2., len(a)-1)
return mean - m*sem, mean + m*sem
Для R я проверил t.test(a). Доверительный интервал GraphPad для средней страницы содержит информацию "пользовательского уровня" о зависимости размера выборки.
Вот вывод для примера Габриэля:
In [2]: a = np.array([1,2,3,4,4,4,5,5,5,5,4,4,4,6,7,8])
In [3]: confIntMean(a, 0.68)
Out[3]: (3.9974214366806184, 4.877578563319382)
In [4]: st.norm.interval(0.68, loc=np.mean(a), scale=st.sem(a))
Out[4]: (4.0120010966037407, 4.8629989033962593)
Обратите внимание, что разница между confIntMean() а также st.norm.interval() интервалы здесь относительно невелики; len(a) == 16 не слишком маленький.
Я проверил ваши методы, используя массив с известным доверительным интервалом. numpy.random.normal(mu,std,size) возвращает массив с центром в mu со стандартным отклонением std (в документах это определяется как Standard deviation (spread or “width”) of the distribution.).
from scipy import stats
import numpy as np
from numpy import random
a = random.normal(0,1,10000)
mean, sigma = np.mean(a), np.std(a)
conf_int_a = stats.norm.interval(0.68, loc=mean, scale=sigma)
conf_int_b = stats.norm.interval(0.68, loc=mean, scale=sigma / np.sqrt(len(a)))
conf_int_a
(-1.0011149125527312, 1.0059797764202412)
conf_int_b
(-0.0076030415111100983, 0.012467905378619625)
Так как значение сигмы должно быть от -1 до 1, / np.sqrt(len(a)) Метод представляется неверным.
редактировать
Поскольку у меня нет репутации комментировать выше, я поясню, как этот ответ связан с полным ответом unutbu. Если вы заполняете случайный массив нормальным распределением, 68% от общего числа попадет в 1-σ от среднего значения. В случае выше, если вы проверите, что вы видите
b = a[np.where((a>-1)&(a <1))]
len(a)
> 6781
или 68% населения попадает в 1σ. Ну около 68%. По мере того, как вы используете все больший и больший массив, вы приближаетесь к 68% (в испытании 10, 9 были между -1 и 1). Это потому, что 1-σ является неотъемлемым распределением данных, и чем больше у вас данных, тем лучше вы сможете их разрешить.
По сути, моя интерпретация вашего вопроса была такова: если у меня есть образец данных, которые я хочу использовать для описания распределения, из которого они взяты, каков метод для определения стандартного отклонения этих данных? в то время как интерпретация УНУТБУ, кажется, больше. Каков интервал, в который я могу поместить среднее с уверенностью 68%?, Что означало бы, для желейных бобов, я ответил, Как они угадывают, и unutbu ответил, Что их догадки говорят нам о желейных бобах.