Правильный способ получения доверительного интервала с помощью scipy
У меня есть 1-мерный массив данных:
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)))
т.е. на сигме используется коэффициент 1/np.sqrt(len(a))
.
Вопрос: какая версия правильная?
Ответы
Ответ 1
68% доверительный интервал для одной ничьей из нормального распределения с среднее значение mu и std sigma
stats.norm.interval(0.68, loc=mu, scale=sigma)
68% -ный доверительный интервал для среднего значения N вытягивает из нормального распределения со средним значением mu и std sigma является
stats.norm.interval(0.68, loc=mu, scale=sigma/sqrt(N))
Интуитивно эти формулы имеют смысл, так как если вы задержите кувшин желе beans и попросите большое количество людей угадать количество желе beans, каждый человек может быть выключен много - такое же отклонение std sigma
, но среднее значение догадок сделает замечательно прекрасную работу по оценке фактического числа, и это отражается стандартным отклонением среднего сокращения в коэффициенте 1/sqrt(N)
.
Если одна ничья имеет дисперсию sigma**2
, то по формуле Bienaymé сумма 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)
в знаменателе.
Вот пример кода, основанного на коде Tom, который демонстрирует приведенные выше утверждения:
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
с желаемым
частота.
Если вы берете образец из дистрибутива и вычисляете
среднее значение образца и отклонение std,
mean, sigma = a.mean(), a.std()
будьте осторожны, чтобы отметить, что нет никакой гарантии, что они
равное среднему значению и стандартным отклонениям, и что мы принимаем
население обычно распределяется - это не автоматические деньги!
Если вы возьмете образец и хотите оценить среднюю и стандартную численность населения
отклонения, вы должны использовать
mean, sigma = a.mean(), a.std(ddof=1)
поскольку это значение для сигмы является несмещенной оценкой для стандартного отклонения населения.
Ответ 2
Я просто проверил, как R и GraphPad вычисляют доверительные интервалы, и они увеличивают интервал в случае небольшого размера выборки (n). Например, более чем в 6 раз для n = 2 по сравнению с большим n. Этот код (на основе shasan answer) соответствует их доверительным интервалам:
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 не слишком мала.
Ответ 3
Я проверил ваши методы, используя массив с известным доверительным интервалом. numpy.random.normal(mu, std, size) возвращает массив с центром в mu со стандартным отклонением std (в docs, это определяется как 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 & sigma; от среднего. В приведенном выше случае, если вы проверите, что видите
b = a[np.where((a>-1)&(a <1))]
len(a)
> 6781
или 68% населения находится в пределах 1 & sigma;. Ну, около 68%. Поскольку вы используете массив большего и большего размера, вы будете приближаться к 68% (в испытании 10, 9 были между -1 и 1). Это потому, что 1- & сигма; является неотъемлемым распределением данных, и чем больше данных у вас есть, тем лучше вы можете его решить.
В принципе, моя интерпретация вашего вопроса была Если у меня есть образец данных, которые я хочу использовать для описания распределения, из которого они были сделаны, то какой метод можно найти для стандартного отклонения этих данных? в то время как интерпретация unutbu представляется более . Каков интервал, на который я могу разместить среднее значение с доверием 68%?. Это будет означать, что для желе beans я ответил. Как они угадывают и unutbu ответил. Что их догадки рассказывают нам о желе beans.