Optim e New-Raphson

Exercício 4)

Gerar 100 valores da distribuição normal(\(\mu\) = 3, \(\sigma ^2\) = 4). Utilizando o método de Newton-Rapson, com cálculo numérico das derivadas, obtenha as estimativas de máxima verossimilhança dos parâmetros \(\mu\) e \(\sigma ^2\). Obtenha também a estimativa adotando a função optim e compare os resultados.

Densidade da Distribuição Normal

\[ f_X(x|\mu,\sigma^2) = \frac{1}{\sqrt{2\pi}\sigma} e^{\frac{-(x-\mu)^2}{2\sigma^2}} \]

Função de Verossimilhança

Achando a verossim. de \(f_X(x)\)

\[ L(\mu,\sigma^2|x) = \prod_{i=1}^{n} f_X(x_i|\mu,\sigma^2) = \prod_{i=1}^{n} \frac{1}{\sqrt{2\pi}\sigma} e^{\frac{-(x-\mu)^2}{2\sigma^2}} \]

\[ L(\mu,\sigma^2|x) = {\left( \frac{1}{\sqrt{2\pi}\sigma} \right)} ^n e^{-\frac{1}{2\sigma^2}\sum_{i=1}^{n}(x_i-\mu)^2} \]

Log-Verossimilhança

Aplicando ln(.) em \(L(\mu,\sigma^2|x)\)

\[ l(\mu,\sigma^2|x) = ln \left( \frac{1}{\sqrt{2\pi}\sigma} \right)^n e^{-\frac{1}{2\sigma^2}\sum_{i=1}^{n}(x_i-\mu)^2} \]

\[ l(\mu,\sigma^2|x) = -nln(\sqrt{2\pi}) - nln(\sigma) - \frac{1}{2\sigma^2}\sum_{i=1}^{n}(x_i-\mu)^2 \]

Derivadas Parciais

Fazendo as derivadas parciais em relação a cada parâmetro

Primeiras derivadas Parciais

Em relação a \(\mu\) :

\[ \frac{\partial l(\mu,\sigma^2|x)}{\partial \mu} = \frac{-n\mu}{\sigma^2} + \frac{\sum_{i=1}^{n}x_i}{\sigma^2} \]

Em relação a \(\sigma ^2\) :

\[ \frac{\partial l(\mu,\sigma^2|x)}{\partial \sigma^2} = \frac{-n}{\sigma}+\frac{\sum_{i=1}^{n}(x_i-\mu)^2}{\sigma^3} \]

O Vetor de Primeiras derivadas parciais fica:

\[g(\mu,\sigma^2) = \begin{pmatrix}\frac{\partial l(\mu,\sigma^2|x)}{\partial \mu} & , & \frac{\partial l(\mu,\sigma^2|x)}{\partial \sigma^2} \end{pmatrix}\]

Segundas Derivadas Parciais

Em relação a \(\mu\):

\[ \frac{\partial^2 l(\mu,\sigma^2|x)}{\partial \mu^2} = \frac{-n}{\sigma^2} \]

Em relação a \(\sigma^2\) é

\[ \frac{\partial^2 l(\mu,\sigma^2|x)}{\partial \sigma^2 \partial \sigma^2} = \frac{n}{\sigma^2} - \frac{3\sum_{i=1}^{n}(x_i-\mu)^2}{\sigma^4} \]

E a mista:

\[ \frac{\partial^2 l(\mu,\sigma^2|x)}{\partial \mu \partial \sigma^2} = \frac{-2(\sum_{i=1}^{n}x_i - n\mu)}{\sigma^3} \]

A Matriz Hessiana fica:

\[H(\mu,\sigma^2) = \begin{pmatrix} \frac{\partial^2 l(\mu,\sigma^2|x)}{\partial \mu^2} & \frac{\partial^2 l(\mu,\sigma^2|x)}{\partial \mu \partial \sigma^2} \\ \frac{\partial^2 l(\mu,\sigma^2|x)}{\partial \sigma^2 \partial \mu} & \frac{\partial^2 l(\mu,\sigma^2|x)}{\partial \sigma^2 \partial \sigma^2} \end{pmatrix}\]

Métodos

Encontrando estimativas de máxima-verossimilhança usando os dois métodos solicitados

Método New-Raphson

Implementando as derivadas parciais no método New-Raphson:

set.seed(2024)
x <- rnorm(100,mean=3,sd=2)
g <- function(x,theta){
  n <- length(x)
  mi<-theta[1]
  sigma2<-theta[2]
  part1 <- (-n*mi)/sigma2 + sum(x)/sigma2
  part2 <- -n/sqrt(sigma2) + (sum((x-mi)^2))/(sigma2*sqrt(sigma2))
  return(c(part1,part2))
}
H <- function(x,theta){
  n<-length(x)
  mi<-theta[1]
  sigma2<-theta[2]
  part1 <- -n/sigma2
  part2 <- (-2*(sum(x)-n*mi))/(sigma2*sqrt(sigma2))
  part3 <- n/(sigma2) - (3*sum((x-mi)^2))/sigma2^2
  h<-matrix(c(part1,part2,part2,part3),ncol=2,nrow=2,byrow=T)
  return(h)
}

n.max<-100
theta<-matrix(NA,nrow=2,ncol=n.max)
dif_<-1
epsilon<-0.001
cont<-1
theta[,1]<-c(1,1)
while(abs(dif_)>=epsilon&cont<=n.max){
  theta[,cont+1]<-theta[,cont]-solve(H(x,theta[,cont]))%*%g(x,theta[,cont])
  cat("theta[",cont+1,"]=",theta[ ,cont+1],"\n")
  dif_<-sqrt(sum((theta[,cont+1]-theta[,cont])^2))
  cont<-cont+1
}

theta[ 2 ]= 2.925056 0.9740673 
theta[ 3 ]= 2.882537 1.246495 
theta[ 4 ]= 2.857256 1.535415 
theta[ 5 ]= 2.843112 1.831865 
theta[ 6 ]= 2.835793 2.126897 
theta[ 7 ]= 2.832348 2.412127 
theta[ 8 ]= 2.830899 2.680279 
theta[ 9 ]= 2.830364 2.925653 
theta[ 10 ]= 2.830194 3.144434 
theta[ 11 ]= 2.830148 3.334774 
theta[ 12 ]= 2.830137 3.496639 
theta[ 13 ]= 2.830135 3.631473 
theta[ 14 ]= 2.830135 3.741756 
theta[ 15 ]= 2.830135 3.830542 
theta[ 16 ]= 2.830135 3.901073 
theta[ 17 ]= 2.830135 3.956489 
theta[ 18 ]= 2.830135 3.999638 
theta[ 19 ]= 2.830135 4.032997 
theta[ 20 ]= 2.830135 4.05864 
theta[ 21 ]= 2.830135 4.078265 
theta[ 22 ]= 2.830135 4.093233 
theta[ 23 ]= 2.830135 4.104619 
theta[ 24 ]= 2.830135 4.113262 
theta[ 25 ]= 2.830135 4.119813 
theta[ 26 ]= 2.830135 4.124772 
theta[ 27 ]= 2.830135 4.128524 
theta[ 28 ]= 2.830135 4.131359 
theta[ 29 ]= 2.830135 4.133501 
theta[ 30 ]= 2.830135 4.135119 
theta[ 31 ]= 2.830135 4.13634 
theta[ 32 ]= 2.830135 4.137262 

Função Optim()

Usando a função optim

set.seed(2024)
x <- rnorm(100,3,2)
f1 <- function(theta){
n <- length(x)
mi <- theta[1]
sigma2 <- theta[2]

 z <- -n*log(sqrt(2*pi*sigma2)) - sum((x-mi)^2)/(2*sigma2)
 return(-z)
 }
 chute <- c(1, 1)
 optim(chute, f1, method="L-BFGS-B",lower=c(-Inf,0),upper=c(Inf,Inf))

$par
[1] 2.830135 4.140099

$value
[1] 212.9298

$counts
function gradient 
      12       12 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"