Разберемся с измерением случайных величин (ч. 4/4)

Доверительный интервал

Доверительный интервал – интервал, который с заданной вероятностью накроет неизвестное значение оцениваемого параметра распределения (вообще говоря, это относится к любому измеримому значению).

Оценка какого-либо параметра распределения с помощью доверительного интервала называют интервальной оценкой.

Главное преимущество интервальной оценки – возможность характеризовать уверенность в вычисленных параметрах распределения. Например, можно будет сказать: «С вероятностью 95% среднее значение \mu нормального распределения для моих данных будет лежать в интервале от 100 до 102, или, как еще можно это записать, \mu=100\pm1.

Зачем вообще нужна интервальная оценка? Предположим, что мы хотим узнать время обработки запроса каким-либо веб-сервисом. Для этого мы производим несколько измерений: мы последовательно посылаем запросы, получаем ответы на запросы, записываем время между отправкой запроса и получением ответа (время отлика). Как вы уже могли заметить из предыдущих записей, в ходе измерения мы получаем множество различных времен отклика. Какие-то значения встречаются чаще, какие-то – реже. Всё это говорит о том, что мы имеем дело со случайной величиной. У случайной величины есть некий закон распределения, который задается рядом параметров. Параметры распределения мы не знаем, но можем их оценить на основе имеющейся выборки (результатов измерений). Например, оценкой дисперсии нормального распределения будет выборочная дисперсия.

Однако, наша оценка – это тоже случайная величина. В этом легко убедиться, если запустить программу сбора данных из первой записи еще раз. Вы получите примерно те же данные, но тем не менее они будут другими, а значит и выборочная дисперсия будет другой. Другими словами эта оценка (или, как еще говорят, статистика) зависит от конкретных выборок.

Согласно центральной предельной теореме при большом количестве измерений наша оценка параметра распределения имеет распределение, близкое к нормальному. Поэтому, чем больше у нас будет выборка, тем ближе выборочная дисперсия будет к дисперсии (параметру) нормального закона распределения! Но когда остановится? Понятно, что если мы сделаем 10 000 измерений, то мы достаточно (с уверенностью в 99,9999%) точно сможем оценить истинные значения параметров распределения, но часто такой возможности нет. На помощь приходит интервальная оценка.

Мы можем сказать так: «Мы не хотим делать 10 000 измерений, но нам достаточно быть уверенным в наших оценках параметров распределения на 95%». Проведя вычисления доверительного интервала, вполне может оказаться, что для 95% уверенности достаточно всего лишь 50 измерений. Существенная экономия сил и ресурсов!)

Обозначим за (1-\alpha)*100\% доверительную вероятность, или надежность, или доверительный уровень (англ. level of confidence) – вероятность того, что истинное значение измеряемой окажется в доверительном интервале. Саму величину \alpha называют уровнем значимости. Математически 1-\alpha можно записать так:
P(a\le \theta \le b)=1-\alpha,
где \theta – это истинное значение оцениваемого интервала, a – нижнее значение интервала, b – верхнее.

Интересно, что \alpha уже встречался ранее, также назывался уровнем значимости, но обозначал вероятность отвержения верной нулевой гипотезы.

Итак, чем больше \alpha, тем меньше доверительная вероятность и тем уже доверительный интервал. Обычно достаточно \alpha=0.05\%, тогда доверительная вероятность будет равна 95%.

У разных распределений разные параметры и, соответственно, построение интервальных оценок для этих параметров разное.

Доверительные интервалы нормального распределения

Посмотрим, как выглядит доверительный интервал для математического ожидания нормального распределения при известном стандартном отклонение \sigma (напоминаю, стандартное отклонение – это квадратный корень от дисперсии). Возникает вопрос, а откуда нам может быть известно стандартное отклонение? Эта характеристика распределения случайной величины может быть задана, например, при производстве какого-либо изделия технический процесс должен обеспечит заданное стандартное отклонение. В этой ситуации доверительный интервал, построенные на базе случайной выборки, поможет определить – соответствует ли изготовленная продукция стандартам тех. процесса.

\bar{X}\pm Z_{1-\alpha}\frac{\sigma}{\sqrt{n}}, где n – объем выборки, \overline{X} – выборочное среднее и, самое интересное, Z_{1-\alpha} – квантиль стандартного нормального распределения (нормальное распределение с мат. ожиданием = 0 и стандартном отклонением = 1) для заданной вероятности 1-\alpha.

Не будем приводить все выкладки, почему именно так определяется доверительный интервал. Подробнее вы можете почитать в [1] (стр. 50-55).

Поясню еще раз про квантиль. Квантиль – это такое значение x, которое случайная величина не превысит с заданной вероятностью 1-\alpha. Квантиль – это значение обратной функции распределения, т.е. 1-\alpha=F(Z).

Статистика малых выборок

В случаях, когда стандартное отклонение неизвестно, а также объем выборки не превышает 10-30 значений, используется t-статистика для оценки значения параметра распределения. t-статистика вычисляется как квантиль t-распределения (распределение Стьюдента, см. первую часть этой серии записей) и задается только одним параметром df – степенью свободы. Для нормального распределения df=n-1, где n – объем выборки.

Тогда интервальная оценка мат. ожидания будет выглядеть следующим образом:

\bar{X}\pm t_{\alpha,df} \frac{S_x^2}{\sqrt{n}}, где t_{\alpha,df} – коэффициента Стьюдента (значение квантили статистики t порядка P=1-\alpha/2 для числа степеней свободы df = n – 1), S_x^2 – выборочная дисперсия.

Доверительные интервалы экспоненциального распределения

После предыдущей записи, я попробовал еще одно распределение для имеющихся данных – экспоненциальное. Оказалось, что оно лучше подходит, чем распределение Эрланга. Результаты измерений вы найдете на иллюстрации ниже и в прилагаемом файле ipython notebook. Видимо, это произошло из-за того, что распределение этих конкретных данных не имеет левого хвоста.

Поэтому, продолжая наш пример, найдем доверительный интервал для параметра \lambda экспоненциального распределения. Кстати, если бы мы использовали экспоненциальное распределение в расчетах в рамках теории массового обслуживания, то параметр \lambda называли «интенсивностью поступления заявок» и он обозначал бы среднее время поступления новых заявок за единицу времени.

Оценкой \lambda является обратное значение среднего арифметического по выборке, т.е.
\hat\lambda=1/\bar{x}.

Доверительный интервал выглядит следующим образом:
$latex \frac{2n}{\widehat{\lambda}\chi^2_{1-\alpha/2,2n}} < \frac{1}{\lambda} < \frac{2n}{\widehat{\lambda}\chi^2_{\alpha/2,2n}}$, где $latex \chi^2_{\alpha/2,2n}$ – квантиль распределения хи-квадрат при заданных доверительной вероятностью и степенью свободы $latex k=2n$.

Код ниже проводит сначала графическую оценку параметра распределения, потом точечную с помощью критерия Пирсона и далее интервальную.

#=========================================
# 1) Графический метод
#=========================================

hist(latencies, bins_number, normed=1)  # Построим нормализованную гистограмму

print(ss.expon.fit(latencies))
exp_loc, exp_scale = ss.expon.fit(latencies)        # Оценим положение распределения

# Наберем больше промежуточных точек для плавности графика
x = np.linspace(np.min(latencies), np.max(latencies) + 10, number_of_samples * 2)

# Получим соответствующие значение по функции плотности вероятности для распределения Эрланга
# В качестве начала построения (параметр loc) выберем минимальное значение  из выборки,
# т.е. сместим распределение так, чтобы оно начиналось с нуля.
# В scipy параметр scale = 1 / lambda, где lambda - параметр экспоненциального распределения.
distr_data = ss.expon.pdf(x, loc=np.min(latencies), scale=4.8)

# Построим график
plot(x, distr_data, 'r--')


#=========================================
# 2) Точечная оценка параметра распределения
#=========================================


# Оценим с помощью хи-критерия соответствие распределения
counts, bins = np.histogram(latencies, bins=bins_number)
p_observed = np.array(counts, dtype=float)

tmp = ss.expon.cdf(bins, loc=exp_loc-0.6, scale=4.8)
#print(tmp)

# Now get probabilities for intervals
p_expected = []
for idx, p in enumerate(tmp):
    if idx == len(tmp) - 1: break
    p_expected.append(tmp[idx+1] - tmp[idx])

p_expected = number_of_samples*np.array(p_expected)

# Calculate using scipy built in chi-square test
chi_test_value, probability_fit = ss.chisquare(p_observed, np.array(p_expected, dtype=float))
print('Вероятность того, что экспоненциальное распределение подходит к полученным данным: {0:f}'.format(probability_fit))

# Calculate by ourselved
chi_star = np.sum((p_observed-p_expected)**2/p_expected)
print("chi_star = {}".format(chi_star))
conf_interval = 0.95
df = len(bins) - 3
chi = ss.chi2.ppf(conf_interval, df)  # chi-square quntile for alpha = conf-interval, degrees of freedom = df
print("chi = {}".format(chi))

Получим следующее:

Следующий код проводит интервальную оценку параметра:

#==================================================
# 3) Интервальная оценка параметра распределения
#==================================================


n = size(latencies)
alpha = 0.05

# Сдвинем данные так, чтобы они начинались с 0
data = np.array(latencies) - np.min(latencies)

# Посчитаем оценку лямбда
lambda_estimate = 1/np.mean(data)
print("Оценка 1/lambda на основе выборки: {}".format(1/lambda_estimate))

# Посчитаем левую границу интервала
a = 2 * n / (lambda_estimate * ss.chi2.ppf(1 - alpha / 2, 2 * n))

# Правую
b = 2 * n / (lambda_estimate * ss.chi2.ppf(alpha / 2, 2 * n))
             
print("95% доверительный интервал: {} < 1/lambda < {}".format(a, b))


ci_width = b - a
print("Ширина доверительного интервала: {}".format(ci_width))

Получим следующее:

    Оценка lambda на основе выборки: 0.249376558603
    Оценка 1/lambda на основе выборки: 4.01
    95% доверительный интервал: 3.32700158323 < 1/lambda < 4.92847012339

Я везде привожу \frac{1}{\lambda}, т.к. это значение использовалось в графическом методе вместо просто \lambda.

Определения необходимого объема выборки для нормального распределения при известном стандартном отклонении

А теперь разберемся с самым интересным! Вот ради чего всё и затевалось.)

До этого мы заранее фиксировали объем выборки, не учитывая доверительную вероятность. В действительности же сначала задается, насколько точной должна быть оценка, и только потом определяют объем выборки.

Опять начнем с простой ситуации, когда известно стандартное отклонение. Такой вариант маловероятен, но тем не менее возможен. Вспомним запись мат. ожидания с доверительным интервалом:
\bar{X}\pm Z_{1-\alpha} \frac{\sigma}{\sqrt{n}}
Величина, добавляемая и вычитаемая из \bar{X}, называется ошибкой выборочного исследования и обозначается следующим образом:
e=Z_{1-\alpha} \frac{\sigma}{\sqrt{n}}

Из формулы e можно определить объем выборки n:
n=\frac{Z^2 \sigma^2}{e^2}

Таким образом, для определения объема выборки необходимо знать три параметра.
1. Требуемый доверительный уровень, который влияет на величину Z, являющуюся критическим значением стандартизованного нормального распределения.
2. Приемлемую ошибку выборочного исследования e. Задается в измеряемых единицах. Например, если бы у нас было нормальное распределение e можно было бы взять равным 5 (мс).
3. Стандартное отклонение ?.

Определения необходимого объема выборки для нормального распределения при неизвестном стандартном отклонении

Здесь, как и при определении доверительного интервала, воспользуемся t-распределением и выборочной дисперсией S_x^2.

n\ge t^2_{\alpha,df}  (\frac{S_x}{\delta})^2=t^2_{\alpha,df} \epsilon^2,
где \delta – требуемая точность. Обычно требуемую точность берут в размерах выборочной дисперсии. Например, если нам нужна точность измерений в 2 раза выше выборочной дисперсии, то \epsilon=\frac{S_x}{0.5S_x}=2.

Стоит отметить, что по сути выборочная дисперсия зависит от объема выборки, поэтому вычисляемый объем выборки необходимо решать методом последовательных приближений. В качестве первого значения выборки можно взять, например, df=10. Посчитать n. Получим, например, 8. Снова n, но при df=8-1=7. Получим n=7. И т.д., пока значение n станет приблизительно неизменным.

Определения необходимого объема выборки для экспоненциального распределения

Вернемся к нашему примеру. Мы установили, что время отлика сайта имеет экспоненциальное распределение. Посчитали доверительный интервал для нашей выборки в 100 значений.

Поставим следующую задачу: сколько необходимо измерений, чтобы быть уверенным в том, что оценка параметра \lambda попадет в доверительный интервал с вероятностью 95%?

Если честно ни русской, ни в иностранной литературе ответа на этот вопрос я не нашел, поэтому решил вывести формулу сам по типу того, как это делали с нормальным распределением.

Из неравенства доверительного интервала выведем формулу расчета необходимого объема выборки. Значение ширины доверительного интервала рассчитываем как разницу между верхней и нижней интервальными оценками и примем за \epsilon (эпсилон). Это значение мы зафиксируем.

\frac{1}{n} = \frac{2}{\epsilon \widehat{\lambda}} (\frac{1}{\chi^2_{\alpha/2,2n}} - \frac{1}{\chi^2_{1-\alpha/2,2n}})

Код ниже просто рассчитывает приведенную формулу:

#==================================================
# 4) Определение необходимого объема выборки 
# при фиксированной ширине доверительного интервала
#==================================================

n = 100
epsilon = ci_width
alpha = 0.1

n_ci = (2 / (epsilon * lambda_estimate)) * (1/ss.chi2.ppf(alpha / 2, 2 * n) - 1/ss.chi2.ppf(1 - alpha / 2, 2 * n))
n_ci = 1 / n_ci
print('Необходимый объем выборки = {} для попадания истинного значения параметра распределения' \
     ' в интервал шириной = {} с доверительной вероятностью = {}%'.format(round(n_ci), epsilon, (1-alpha)*100))

Результ: «Необходимый объем выборки = 120.0 для попадания истинного значения параметра в интервал шириной = 1.60146854016 с доверительной вероятностью = 90.0%».

Нахождение объема выборки при фиксированной доверительной вероятности и ширине доверительного интервала

Дальше идут весьма сомнительные рассуждения, поэтому не рекомендую их читать (тем более всерьез). =) Задача в заголовке раздела вполне может появится, но весьма редко. В ходе решения предыдущей задачи (когда мы могли изменять доверительную вероятность) я неспециально решил задачу в заголовке.

Как в случае с нормальным распределением, формула вышла рекуррентной. Оптимальное значение можно найти с помощью метода последовательных приближений или с помощью линейной регрессии. Для использования метода линейной регрессии необходимо вычислить частную производную по n, что в нашем случае не так-то просто, поэтому ограничимся простым методом приближения.

Код (который вы также найдете в файле ipython notebook в конце поста):

#==================================================
# 5*) Определение необходимого объема выборки 
# при фиксированной ширине доверительного интервала и доверительной вероятности
#==================================================
alpha = 0.05
epsilon =  0.5 * ci_width
optimum_n = 0
N = [0]

for i in range(0, 500, 10):
    n = i
    
    for k in range(5):
        n_ci = (2 / (epsilon * lambda_estimate)) * (1/ss.chi2.ppf(alpha / 2, 2 * n) - 1/ss.chi2.ppf(1 - alpha / 2, 2 * n))
        n_ci = 1 / n_ci
        n = n_ci
        N.append(n_ci)
        #print("{}".format(n_ci))
    
    #print('---')
    
    # Замечено, что при увеличении первоначального значения n
    # в определенной итерации метода приближения вычисленный объем выборки
    # перестает падать или расти, а колеблется около одного и того же числа.
    # Это число можно считать оптимальным объемом выборки при заданной ширине доверительного интервала.
    # Что требует доказательства.
    
    # Ниже мы грубо оцениваем размах колебания значения n_ci в этой итерации.
    # Если оно мало, значит это значение n можно считать оптимальным.
    if np.std(N[-5:]) < 5:
        optimum_n = n_ci

print('Эверстически оптимальный объем выборки: {}'.format(math.ceil(optimum_n)))
print('При заданной ширине доверительного интервала: {}'.format(epsilon))
print('И доверительной вероятности: {}%'.format((1-alpha)*100))
plot(N)

Если использовать в качестве начального значения n изначальный объем выборки (в нашем случае это 100), то каждое следующее приближение будет уменьшать вычисленной значение n_ci до нуля, однако, я заметил, что при увеличении первоначального значения n (больше первоначального объем) в определенной итерации метода приближения вычисленный объем выборки перестает падать или расти, а колеблется около одного и того же числа. Это число можно считать оптимальным объемом выборки при заданной ширине доверительного интервала. Это, конечно, требует доказательства и является полнейшей эвристикой.

Полученный с помощью этой эвристики результаты хорошо согласуются с предполагаемым поведением этой функции: чем шире значение доверительного интервала, тем меньше нужен объем выборки и наоборот. Однако сам алгоритм нестабилен и не выдает результатов при \epsilon на порядок меньше исходной ширины интервала. При больших значениях он соответственнно выдаст 0 (что можно понимать как выборку с 1 значением), т.к. если интервал достаточно широк, то попасть в него будет очень просто с любой выборкой.

Ниже график, на котором показано вычисление объема выборки методом последовательных приближений при различных начальных значениях n. Хорошо виден интервал, в котором вычисленное значение колеблется, но заметно не растет и не падает.

В общем, если что-то подобное пригодится в будущем, я спрошу совета у продвинутых математиков.) Если тот, кто читает эту запись, предложит другой способ определения объема выборки для экспоненциального распределения, я с удовольствием его послушаю!

Скачать файл с ipython notebook раздела

Заключение

Занимаясь написанием этой череды постов, я узнал много нового про распределения и методы работы с ними, надеюсь и вы почерпнули для себя что-нибудь интересное и полезное. =)

Список литературы

  1. Спирин H.A. Методы планирования и обработки результатов инженерного эксперимента : Конспект лекций. / H.A. Спирин, В.В. Лавров. Екатеринбург : ГОУ ВПО УГТУ-УПИ, 2004. — 257 с.
  2. Статистика для менеджеров с использованием Microsoft Office Excel 4-e издание. Дэвид М. Левин и др. – отлично рассказано в 7 главе про практическое использование доверительных интервалов. И именно 7-я глава есть в открытом доступе.
  3. Confidence Intervals, Guy Lebanon – отличная работа, рассказывающая о доверительных интервалах.
  4. http://en.wikipedia.org/wiki/Exponential_distribution#Maximum_likelihood – про оценку параметра экспоненциального распределения.