Данный ноутбук является конспектом по курсу «R для теории вероятностей и математической статистики» (РАНХиГС, 2017-2018). Автор ноутбука вот этот парень по имени Филипп. Если у вас для него есть деньги, слава или женщины, он от этого всего не откажется. Ноутбук распространяется на условиях лицензии Creative Commons Attribution-Share Alike 4.0. При использовании обязательно упоминание автора курса и аффилиации. При наличии технической возможности необходимо также указать активную гиперссылку на страницу курса. На ней можно найти другие материалы. Фрагменты кода, включенные в этот notebook, публикуются как общественное достояние.
Метод максимального правдоподобия — это основная лошадка современной статистики. Функция правдоподобия встречается абсолютно во всех областях от эконометрики и байесовской статистики до нейронных сетей и рекомендательных систем. Из-за этого очень важно как следует понять как она устроена.
В этой тетрадке именно этим мы и займёмся. Для начала мы попробуем посмотреть на поднаготную метода максимального правдоподобия и понять где это там в нём зарыта информация Фишера. И почему это вообще информация. Скорее всего, вас довольно давно будоражит этот вопрос. Затем мы решим в R парочку задачек.
library("ggplot2") # Пакет для красивых графиков
library("grid") # Пакет для субплотов
# Отрегулируем размер картинок, которые будут выдаваться в нашей тетрадке
library('repr')
options(repr.plot.width=4, repr.plot.height=3)
Пришло время формул! Итак, пусть в нашем распоряжении оказалась выборка $x_1, \ldots x_n$. Пусть мы предположили, что эту выбору на нас из сундука выплюнула случайная величина с плотностью распределения $f(x \mid \theta)$. Про параметр $\theta$ нам ничего не известно и мы хотели бы оценить его по выборке.
Запишем правдоподобие наших данных, то есть вероятность пронаблюдать именно эту выборку:
$$L( \theta \mid x_1, \ldots, x_n) = f(x_1, \ldots, x_n \mid \theta) = \prod_{i=1}^n f(x_i \mid \theta)$$Наша выборка - фиксирована. При разных значениях $\theta$ мы получаем большую или меньшую вероятность получить выборку близку к наблюдаемой. Было бы круто найти такое значение $\theta$, для которого эта вероятность оказалась бы максимальной.
Обратите внимание, что функцию $L( \theta \mid x_1, \ldots, x_n)$ уже нельзя трактовать как плотность распределения. Параметр $\theta$ здесь не рассматривается как случайная величина, он фиксирован. Интеграл $L( \theta \mid x_1, \ldots, x_n)$ по всем возможным значениям $\theta$ не будет равен единице.
Для максимизации правдоподобия нам придётся брать производные. Куда приятнее их брать от логарифма, потому что все произведения превратятся в суммы и жить станет проще. Выпишем логарифм правдоподобия:
$$ \ln L( \theta \mid x_1, \ldots, x_n) = \sum_{i=1}^n \ln f(x_i \mid \theta). $$Берём производную и приравниваем её к нулю.
$$ \frac{\partial \ln L}{\partial \theta} = \sum_{i=1}^n \frac{\partial \ln f(x_i \mid \theta)}{\partial \theta} = 0 $$Решив это уравнение мы получим оценку максимального правдоподобия.
Посмотрим на одно слагаемое $\ln f(x_i \mid \theta)$. Его можно интерпретировать как логарифм правдоподобия, посчитанного на основании только одного наблюдения $x_i$.
Когда у нас в выборке появляются новые наблюдения, к нашей функции плюсуется дополнительное слагаемое. Интересно было бы посмотреть на то, как при этом функция изменяется.
У нас есть R. Сгенерируем выборку объёма $5$ из нормального стандартного распределения и посмотрим.
x <- rnorm(5)
x
Выпишем для нормального распределения логарифмическую функцию правдоподобия. Она будет зависеть от двух векторов: выборки и параметров.
$$ \ln L(\mu, \sigma^2 \mid x) = -\frac{n}{2}\ln(2\pi)-\frac{n}{2}\ln(\sigma^2)-\sum_{i=1}^n\frac{(x_i - \mu)^2}{2\sigma^2} $$Предположим, что $\sigma$ известна и нам хочется оценить только $\mu$. Перебьём такую функцию в R.
# mu - оцениваемый параметр
# data - вектор данных
log_lik <- function(mu,data) {
n <- length(data)
result <- -n/2*log(2*pi)-n/2*log(1)-sum((data-mu)^2)/(2*1)
return(result)
}
Найдём значение нашей функции в какой-нибудь точке чисто ради фана.
log_lik(1,x)
log_lik(10,x)
Посмотрим как функция выглядит на картинке.
mu = seq(-12,11,0.01)
lnL = rep(0,length(mu))
for(i in 1:length(lnL)){
lnL[i] = log_lik(mu[i],x)
}
qplot(mu, lnL, geom='line')
Выглядит прелестно. Её максимум достигается в районе нуля, но не чётко в нуле. Построим такие же функции для отдельных наблюдений.
mu = seq(-12,11,0.01)
lnL = rep(0,length(mu))
lnL1 = rep(0,length(mu))
lnL2 = rep(0,length(mu))
lnL3 = rep(0,length(mu))
lnL4 = rep(0,length(mu))
lnL5 = rep(0,length(mu))
for(i in 1:length(lnL)){
lnL[i] = log_lik(mu[i],x)
lnL1[i] = log_lik(mu[i],x[1])
lnL2[i] = log_lik(mu[i],x[2])
lnL3[i] = log_lik(mu[i],x[3])
lnL4[i] = log_lik(mu[i],x[4])
lnL5[i] = log_lik(mu[i],x[5])
}
df_L = data.frame(mu=mu,lnL = lnL, lnL1 = lnL1, lnL2 = lnL2,
lnL3 = lnL3,lnL4 = lnL4,lnL5 = lnL5)
ggplot(df_L, aes(x=mu))+
geom_line(aes(y=lnL))+
geom_line(aes(y=lnL1),color='red')+
geom_line(aes(y=lnL2),color='red')+
geom_line(aes(y=lnL3),color='red')+
geom_line(aes(y=lnL4),color='red')+
geom_line(aes(y=lnL5),color='red')
Черная функция для всей выборки равна сумме логарифмических правдоподобий для отдельных наблюдений (красные линии). Она имеет более выраженный максимум по сравнению с логарифмическими функциями для отдельных наблюдений.
Можно построить аналогичную картинку, на которой мы будем постепенно накапливать наблюдения внутри нашей суммы. При добавлении всё новых слагаемых, максимум будет становиться более чётким.
mu = seq(-12,11,0.01)
lnL = rep(0,length(mu))
lnL1 = rep(0,length(mu))
lnL2 = rep(0,length(mu))
lnL3 = rep(0,length(mu))
lnL4 = rep(0,length(mu))
lnL5 = rep(0,length(mu))
for(i in 1:length(lnL)){
lnL[i] = log_lik(mu[i],x)
lnL1[i] = log_lik(mu[i],x[1])
lnL2[i] = log_lik(mu[i],x[1:2])
lnL3[i] = log_lik(mu[i],x[1:3])
lnL4[i] = log_lik(mu[i],x[1:4])
}
df_L = data.frame(mu=mu,lnL = lnL, lnL1 = lnL1, lnL2 = lnL2,
lnL3 = lnL3,lnL4 = lnL4)
ggplot(df_L, aes(x=mu))+
geom_line(aes(y=lnL))+
geom_line(aes(y=lnL1),color='red')+
geom_line(aes(y=lnL2),color='red')+
geom_line(aes(y=lnL3),color='red')+
geom_line(aes(y=lnL4),color='red')
Итак, чем более выпукла функция, тем более выражен максимум, тем более сильно мы уверенны в том, что оценка параметра была найдена хорошо.
Что отвечает за выпуклость функции? Правильно! Вторая производная. Именно её, взятую со знаком минус, интерпретируют как наблюдённую информацию.
$$ I_o(\theta) = - \left( \frac{\partial^2 \ln L}{\partial \theta^2} \right) $$Почему со знаком минус? Потому что у нас максимум, в точке максимума вторая производная отрицательна. А информация должна быть положительной. Если параметр векторный, то в этом случае наблюдэнная информация представляется наблюдённой информационной матрицей:
$$ I_o(\theta) = - \left( \frac{\partial^2 \ln L}{\partial \theta_i \partial \theta_j} \right) = - H $$Индекс o в данном случае обозначает observed, информацию, которую мы реально видели. В матрице на во всех клетках стоят числа. Они зависят от конкретных значений наблюдений.
Обратите внимание, что на первой картинке, правдоподобия, построенные для разных наблюдений выпуклы по разному. Получается, что каждое наблюдение вносит в заострённость нашего правдоподобия разный вклад, то есть в каждом наблюдении содержится разное количество информации.
Математическое ожидание этой матрицы по распределению $x$ называется информационной матрицей Фишера.
$$ I(\theta) = - E\left( \frac{\partial^2 \ln L}{\partial \theta_i \partial \theta_j} \right) = -E(H) $$Ожидаемая информация зависит только от закона распределения наблюдений. Она отражает то, какую информацию в среднем вносит в наше правдоподобие, каждое наблюдение.
Если функция плотности $f(x \mid \theta)$ удовлетворяет условиям регулярности, то тогда для любой несмещённой оценки $\hat \theta$ выполняется неравенство Рао-Фреше-Крамера:
$$ Var(\hat \theta) \ge [I(\theta)]^{-1} $$Более того, в этом случае
$$ I(\theta) = - E\left( \frac{\partial^2 \ln L}{\partial \theta^2} \right) = E\left[ \left( \frac{\partial \ln L}{\partial \theta} \right)^2 \right], $$но доказывать это здесь, мы конечно же не будем.
При тех же условиях регулярности, оценка максимального правдоподоия обладает набором приятных асимптотических свойств:
Асимптотическую нормальность обычно используют на практике для строительства доверительных интервалов и проверки гипотез. Чтобы сделать это, yужно найти подходящую оценку для информации Фишера, посчитанную на основе выборки.
Обычно поступают двумя путями: либо заменяют в матрице Гессе, $-H$, все $\theta$ на оценки, либо заменяют в матрице $-E(H)$ все $\theta$ на оценки. Именно так мы с вами поступим, когда будем решать задачки.