3

estou simulando valores da distribuição Birnbaum-Saunders com parâmetros alpha e beta e posteriormente tento estimar estes valores usando a função optim, apesar da simulação convergir, estou recebendo o erro abaixo, como proceder para corrigir este erro:

n<-1000
z<-rnorm(n,0,1)
alpha=2
beta=2
t<-beta*((alpha*z/2)+sqrt((alpha*z/2)^2+1))^2
remove(par,alpha,beta)  

Loglik<-function(par){
  alpha=par[1]
  beta=par[2]
  ll<-sum(log(t+beta))-n*log(2*alpha)-(n/2)*(log(2*pi*beta))-(3/2)*sum(log(t))-((1/(2*alpha^2))*sum((t/beta)+(beta/t)-2))
  return(-ll)
}
alpha_0=1
beta_0=5
start=c(alpha_0,beta_0)
optim(start,fn=Loglik,method="BFGS",hessian=T)$par

Resultado:

[1] 2.018785 2.133996
There were 25 warnings (use warnings() to see them)

Os erros são:

> warnings()
Mensagens de aviso:
1: In log(t + beta) : NaNs produced
2: In log(2 * pi * beta) : NaNs produced
3: In log(t + beta) : NaNs produced

Já conferi minha verossimilhança, que está correta, não encontro o erro no código?

2 Respostas 2

5

Isso não são erros propriamente ditos, são avisos (warnings) de que algo ocorreu de uma forma que provavelmente não deveria, mas que o código foi executado completamente.

No seu caso, eles são dados porque a função tentou tirar o logaritmo de um número negativo. Pra sua sorte, digamos assim, a função optim consegue lidar bem com isso e ignora o problema, tanto que o resultado obtido parece correto.

De qualquer forma, você pode fazer algumas coisas para evitar os warnings:

  1. Você pode usar supressWarnings() para esconder os avisos. Eles continuarão lá, mas ocultos:

suppressWarnings(optim(start,fn=Loglik,method="BFGS",hessian=T)$par)
  1. Você pode evitar que os valores negativos sejam utilizados no cálculo de ll, adicionando a seguinte linha antes da operação:

if (any(c(t, beta, alpha, pi) < 0)) return(NA)

Note que essas "soluções" são apenas maquiagens para evitar os warnings. Uma solução mais elegante (e provavelmente mais correta) involveria analisar qual motivo de o algoritmo estar gerando valores de alpha e beta negativo, e decidir se isso por si só é um problema (parece que não), para fazer uma escolha mais consciente, ou então resolver o problema na fonte (no método de otimização ou na função sendo otimizada).

1
  • @fsblajinha Que bom, a solução matemática parece a melhor saída. Você pode colocar como resposta, mesmo a pergunta sendo sua!
    – Molx
    Commented 21/05/2016 às 1:24
2

Gostaria de ressaltar outra forma que aprendi hoje que pode também ajudar outras pessoas que é uma dica interessante.

Substituir dentro da função log de verossimilhança os parâmetros alpha e beta por exp(lalpha) e exp(lbeta), e posteriormente, fazer lalpha_0=log(alpha_0) e lbeta_0=log(beta_0). O resultado será exp(par[1]) e exp(par[2]).

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .