Подтвердить что ты не робот

Правильный способ получения доверительного интервала с помощью 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)).

Вопрос: какая версия правильная?

4b9b3361

Ответ 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.