Estimação de Máximo Verossimilhança

Objetivo

Obter estimativas de máximo verossimilhança dos parâmetros da distribuição Birnbaum-Saunders.

Geração de Amostra da Birnbaum-Saunders

  1. Geramos uma amostra de tamanho \(n\) da normal padrão;

  2. Substituímos na fórmula 1) abaixo, gerando uma amostra de tamanho n da distribuição Birnbaum-Saunders.

\[ T=\beta\bigg[\frac{\alpha}{2}Z + \sqrt{\bigg[\frac{\alpha}{2}Z\bigg]^2+1}\;\;\bigg]^2 \]

#=========================== 1o Modo================================================#

alpha<-0.5 ; beta<-2 # valores que o prof. pediu para considerar.

n<-100  # escolha a qtd de valores a ser gerado (isto é, tamanho da amostra)

z<-rnorm(n,mean=0,sd=1)  # n valores da normal padrao

t<-beta*( (alpha/2)*z + sqrt( ((alpha/2)*z)^2 + 1 ) )^2

#============================ 2o Modo===============================================#

amostra.T<-function(n,alpha,beta){
            
                z<-rnorm(n,mean=0,sd=1)  # n valores da normal padrao
  
                t<-beta*( (alpha/2)*z + sqrt( ((alpha/2)*z)^2 + 1 ) )^2
          
          return(t)
}

# t<-amostra.T(100,0.5,2)

Construção da Matriz Hessiana e do Vetor Score ou Gradiente

A matriz Hessiana e o vetor Score são, respectivamente:

\[ H=\begin{pmatrix} H.aa & H.ab \\ H.ba & H.bb \end{pmatrix} \qquad \mbox{e}\qquad U=(U.a\;,\;U.b) \]

em que, \(U.a\), \(U.b\), \(H.aa\), \(H.ab\), \(H.ba\) e \(H.bb\) foram dados.

### Funcao Matriz

Matriz<-function(alpha,beta){
             
             n<-100 
  
             U.a<- -2*n/alpha^3 + alpha^(-3) * sum( t/beta + beta/t ) - n/alpha
  
             U.b<- ( 1/(2*alpha^2) ) * sum( t/beta^2 - 1/t ) - n/(2*beta) + sum( 1/(t+beta) )
  
             Grad<-c(U.a, U.b)
  
               
             H.aa <- (6*n/alpha^4) - (3/alpha^4)*sum( (t/beta) + (beta/t) ) + n/alpha^2 
  
             H.ab<-H.ba<- alpha^(-3) * sum( 1/t - t/beta^2) 
  
             H.bb <-  ( -1/(alpha^2*beta^3) )*sum(t) + n/(2*beta^2) - sum( 1/(t+beta)^2 )  
  
             Hess<-matrix( c(H.aa, H.ab, H.ba, H.bb), nrow = 2, byrow = T )
  
return(list(grad=Grad, hess=Hess))
}

Estimação via Método de Newton

\[\mathbf{\theta^{(t+1)}=\theta^{(t)}-H^{-1}\big(\theta^{(t)}\big)\;U\big(\theta^{(t)}\big)}. \]

alpha<-0.5 ; beta<-2 # parametros a serem estimados

### Estimação ( Método de Newton)

estimacao<-function(alpha0,beta0,tol) {
   
      # valores inicias
          
      theta0<-c(alpha0,beta0) 

      theta<-numeric(2) ; theta<-theta0 ; N<-0
                         
      # Eq. de atualizacao e criterio de parada      
                     
      while((abs(alpha-theta[1])+abs(beta-theta[2])/(abs(theta[1])+abs(theta[2]))) > tol){
                
                Mat<-Matriz(theta[1],theta[2])
    
                theta <- theta - (1/2)*solve( Mat$hess ) %*%  Mat$grad   
                         
                theta<-theta
                
                N<-N+1
            }
  return( setNames(c(theta[1],theta[2],N),c("est. alpha"," est. beta"," No. Interações") ))
}

Exemplo

Neste exemplo, consideramos os valores iniciais \(\alpha_0=0.1\) e \(\beta_0=1\), com tolerância \(\epsilon=0.06\).

Desejamos estimar os parâmetros \(\alpha=0.5\) e \(\beta=2\).

# parametros a serem estimados alpha=0.5 e beta=2 ( os utilizados para gerarem  a amostra)
estimacao(0.1,1.0,0.06)
##      est. alpha       est. beta  No. Interações 
##       0.4642377       2.0584236      15.0000000