Метод максималне веродостојности представља један од статистичких метода за налажење оцена непознатих праметара. Овом методом се добијају оцене које имају оптималне особине. Прецизније, добијају се довољне оцене увак кад оне постоје, као и непристрасне оцене са, бар асимптотски, најмањом дисперзијом.

Нека је \((X_1, X_2, \ldots, X_n)\) случајан узорак из фамилије допустивих расподела \(\{f(x, \theta), \theta \in \Theta\}\) за обележје \(X\).

Фукнкција веродостојности се дефинише као функција од \(\theta \in \Theta\)

\[L(\theta)=L(X_1, X_2, \ldots, X_n;\theta)=\prod_{k=1}^nf(X_k, \theta)\].

За произвољан реализовани узорак \((x_1, x_2, \ldots, x_n)\) при фиксираној вредности параметра \(\theta\), функција

\[L(\theta)=L(x_1, x_2, \ldots, x_n;\theta)=\prod_{k=1}^nf(x_k, \theta)\]

представља вероватноћу реализације узорка \((x_1, x_2, \ldots, x_n)\) и одређује расподелу случајног узорка \((X_1, X_2, \ldots, X_n)\).

Оцена добијена методом максималне веродостојности је она вредност параметра \(\theta=t(x_1, x_2, \ldots, x_n)\) за коју функција веродостојности достиже супремум.

Ако је функција веродостојности узорка диференцијабилна по \(\theta\), оцењена вредност параметра \(\theta\) добијена методом максималне веродостоности одерђује се из једначине веродостојности

\[\frac{\partial}{\partial \theta}L(x_1, x_2, \ldots, x_n;\theta)=0.\]

За вредност \(\theta=t(x_1, x_2, \ldots, x_n)\) која се добија као решење једначине веродостојности треба проверити да ли је то заиста тачка максимума функције веродостојности (проверити да ли је \(\frac{\partial^2}{\partial \theta^2}L\ |_{\theta=t}<0)\).

Како функције \(L(\theta)\) и \(\ln L(\theta)\) због монотоности достижу максимум за исту вредност параметра \(\theta\), често је лакше ову вредност наћи као решење једначине

\[\frac{\partial}{\partial \theta}\ln L(x_1, x_2, \ldots, x_n;\theta)=0.\]

Међутим, у многим ситуациjама jе тешко наћи аналитички облик за оцену максималне веродостоjности. У тим ситуациjама, обично jе потребно израчунати оцену максималне веродостоjности нумеричким путем.

Неке од нумеричких метода које се могу применити су:

\(\color{lightseagreen}{\text{Њутн-Рафсонов алгоритам}}\)

Њутн-Рафсонов алгоритам се користи за налажење решења нелинеарних једначина или система нелинеарних једначина. Видели смо да приликом одређивања оцене методом маскималне веродостојности наилазимо на баш такву једначину, а у случају оцењивања више непознатих параметара треба решити систем.

Претпоставимо да треба наћи решење једначине \(g(x_0)=0\), где је \(g\) диференцијабилна функција. За број \(x\) из околине тачке \(x_0\), применом Тејлоровог развоја око \(x\) се добија \[0=g(x_0)\approx g(x_0) + g^{'}(x)(x_0-x),\] одакле следи \[x_0 \approx x-\frac{g(x)}{g^{'}(x)}.\] Следи да, за дату оцену \(x_k\), можемо одредити оцену \(x_{k+1}\) \[x_{k+1} = x_k -\frac{g(x_k)}{g^{'}(x_k)}.\]

Оваква процедура може да се понови за \(k=1, 2, 3, \ldots\) док \(|g(x_k)|\) (или \(\frac{g(x_k)}{g^{'}(x_k)}\)) не постане довољно мало.

Претпоставимо да је \(\widehat\theta\) оцена добијена методом максималне веродостојности која задовољава једначину \(S(\widehat\theta)=0\), где је \(S(\theta)\) први извод логаритмоване функције веродостојности.

Нека је \(\widehat \theta^{(k)}\) оцена параметра \(\theta\) добијена после \(k\) итерација алгоритма. Тада \[\widehat \theta^{(k+1)}=\widehat \theta^{(k)}+\frac{S(\widehat \theta^{(k)})}{H(\widehat \theta^{(k)})},\] где је \[H(\theta)=-\frac{\partial^2}{\partial \theta^2}\ln L(\theta).\]

Поступак се понавља док \(|S(\widehat \theta^{(k)})|\) или разлика \(|\widehat\theta^{(k)} − \widehat\theta^{(k+1)}|\) не постану довољно мали.

За Њутн-Рафсонов алгоритам је потребно задати почетну (инициjалну) оцену \(\widehat \theta^{(0)}\).

Уколико нисмо сигурни да ли је решење једначине \(S(\widehat\theta)=0\) заиста максимум, иницијалну оцену \(\widehat \theta^{(0)}\) можемо добити методом момената или неком другом методом.

\(\color{lightseagreen}{\text{Пример}}\)

Генеришимо случајан узорак из Кошијеве расподеле са параметром 10.

set.seed(1)
x<-rcauchy(100,location=10)
hist(x) 

median(x)
## [1] 10.13036

Видимо да се параметар локације може приближно оценити медијаном.

У наставку ћемо показати како се на основу Њутн-Рафсоновог метода може оценити параметар локације на основу датог узорка (како у овом случају знамо да је узорак генерисан из Кошијеве расподеле чији је параметар 10, очекујемо да вредност буде приближно 10).

Наћи ћемо нулу првог извода логаритмоване функције веродостојности, која се иначе не може експлицитно изразити.

Нацртајмо функцију чији максимум тражимо.

cauchy_loglike<-function(theta,x)
{
  return(sum(-log(1+(x-theta)^2))-length(x)*pi)
}
n<-50
loglike_array<-vector()
theta<-seq(from=8,to=12,length.out=50)
for(i in 1:n)
{
  loglike_array<-c(loglike_array,cauchy_loglike(theta[i],x))
}

plot(theta,loglike_array,type="l")

Цртање ове функције може бити корисно за избор иницијалне вредности алгоритма.

Направићемо функције које ће рачунати вредности функција \(S\) и \(H\) у итерацијама.

cauchy_score<-function(theta,x)
{
  return(sum(2*(x-theta)/(1+(x-theta)^2)))
}


d_cauchy_score<-function(theta,x)
{
  return(2*sum(((x-theta)^2-1)/(1+(x-theta)^2)^2))
}

Можемо нацртати фукнцију \(S\) и на основу тога видети где је отприлике нула извода, али то може бити тачка локалног максимума или минимума логаритмоване функције веродостојности.

У нашем примеру, пошто параметар можемо лако узорачки оценити - медијаном, то ћемо користити као почетну вредност.

eps<-0.003 # задајемо тачност
theta<-vector() 
theta0<-median(x) #почетна вредност
theta[1]<-theta0

i=1
repeat
{
  theta[i+1]<-theta[i]- cauchy_score(theta[i],x)/d_cauchy_score(theta[i],x)
  i=i+1
  if(abs(theta[i]-theta[i-1])<eps) break
}

Коначно, оцена параметра локације Кошијеве расподеле добојена нумерички методом максималне веродостојности је:

theta[i]
## [1] 9.90179
\(\color{lightseagreen}{\text{Фишеров алгоритам}}\)

Фишеров алгоритам се добија једноставном модификацијом Њутн-Рафсоновог алгоритма. Нека је \({\bf \theta{}}= (\theta_1, \ldots,\theta_p)^T\). У овом алгоритму се \(H\) мења са \(H^*\), где је елемент \((i,j)\) од \(H^*\) дат са \[H^*_{ij}({\bf \theta})=E_\theta[H_{ij}({\bf \theta})]= -E_\theta \left[\frac{\partial^2}{\partial \theta_i\partial \theta_j}\ln f(\bf{X}; {\bf \theta})\right]\] и представља очекивану вредност под претпоставком да је \(\bf{\theta}\) права вредност параметра.

\(H\) се назива Фишерова информациона матрица, а \(H^{*}\) очекивана Фишерова информациона матрица.

Ако jе \(\widehat{\bf \theta}^{(k)}\) оцена параметра \(\bf{\theta}\) након \(k\) итерациjа, тада важи \[\bf{\widehat \theta}^{(k+1)}=\bf{\widehat \theta}^{(k)}+[H(\bf{\widehat \theta}^{(k)})]^{-1}S(\bf{\widehat \theta}^{(k)}).\]

\(\color{lightseagreen}{\text{Пример}}\)

Нека су \((X_1, \ldots, X_n)\) независне jеднако расподељене случаjне величине са Кошиjевом расподелом. Очекивана вредност Фишерове информационе функције је \(\frac{n}{2}\) (израчунати!).

eps<-0.0001 #тачност
theta<-vector() 
theta0<-median(x)
theta[1]<-theta0

i=1
repeat
{
  theta[i+1]<-theta[i]+ cauchy_score(theta[i],x)*2/n
  i=i+1
  if(abs(theta[i]-theta[i-1])<eps) break
}

Добијена је оцена:

theta[i]
## [1] 9.901808