Error in output of Optim when creating a regression structure

Question

By running the simulation below I get the error at the end of the code, someone can help me to solve it?

#Valores iniciais dos parâmetros usados para gerar t
#Ou seja, Valor verdadeiro dos parâmetros
beta0=2
beta1=1
alpha0=3
alpha1=1
truevalue=c(beta0,beta1,alpha0,alpha1)
#Tamanho das amostras
n=100
#Vetor de parâmetros
beta=matrix(c(beta0,beta1),nrow=2,ncol=1)
alpha=matrix(c(alpha0,alpha1),nrow=2,ncol=1)
#Vetor de 1's
const1 <- rep(1,n);
const <- cbind(const1);
#Vetor de cováriavel com distribuição unif(0,1)
X1=matrix(runif(n, 0, 1),nrow=n,ncol=1)
Z1=matrix(runif(n, 0, 1),nrow=n,ncol=1)
#Matrizes de covariáveis
X <-matrix(c(const,X1),nrow=n,ncol=ncol(X1)+1)
Z <-matrix(c(const,Z1),nrow=n,ncol=ncol(Z1)+1)
#Número de colunas de X e Z
p=ncol(X)
q=ncol(Z)
#Vetor de médias
mu=exp(X%*%beta)
#Vetor de dispersão
phi=exp(Z%*%alpha)

#Gerando uma variável aleatória t com distribuição Birnbaum-Saunders(mu,phi)
remove(beta,beta0,beta1,alpha,alpha0,alpha1)

z<-cbind(rnorm(n,0,1))
t<-((phi*mu)/(phi+1))*((z/sqrt(2*phi))+(sqrt((z/sqrt(2*phi))^2+1)))^2

Loglik<-function(theta,dados){
  p=ncol(X)
  q=ncol(Z)
  beta=theta[1:p]
  alpha=theta[p+1:p+q]
  mu=exp(X%*%beta)
  phi=exp(Z%*%alpha)
  lv=sum((phi/2)-(3/2)*log(t)+(1/2)*log(phi+1)-log(4*sqrt(pi*mu))+log(t+(phi*mu/(phi+1)))-(phi/4)*((t*(phi+1)/(phi*mu))+(phi*mu/(t*(phi+1)))))
  #lv=sum(log(((exp(phi/2)*sqrt(phi+1))/(4*sqrt(pi*mu)*t^(3/2)))*(t+((phi*mu)/(phi+1)))*exp((-phi/4)*((t*(phi+1)/(phi*mu))+(phi*mu/(t*(phi+1)))))))
  return(-lv)
}
#Chute inicial para as funções de estimação
start=cbind(1,2,2,3)
#Estimation with function optim
bs_op=optim(start,Loglik,method="BFGS",dados=t,hessian = T) 
bs_op$par

Error in Optim(start, Loglik, method = "BFGS", data = t, Hessian = T): initial value in 'vmmin' is not Finite

Answer 1

It seems that Eduardo’s suggestion in the comments has to do with your problem, yes. In that reply for that question, the hint was given to check the result of its function with the initial values. In your case:

> Loglik(start, t)
[1] NA

This indicates that something is wrong in its function. By investigating it line by line, we can observe a problem here:

> alpha=theta[p+1:p+q]
> alpha
[1] NA NA

The NAonly appear because the account p+1:p+q results in the vector (5, 6), whereas theta only has 4 elements. On second thought, my guess is that you got the operator priority wrong :, and that’s why you got this result. I imagine you wanted the values c(3, 4), and in this case should change the code to:

alpha <- theta[(p+1):(p+q)]

See that now the result is not NAs:

> alpha
[1] 2 3

And the function optim wheel perfectly:

> bs_op <- optim(start,Loglik,method="BFGS",dados=t,hessian = T)
> bs_op$par
        [,1]     [,2]     [,3]     [,4]
[1,] 1.00167 2.078868 1.655399 3.696152

I leave to you the evaluation of this result, because you know better what to expect from so much math. :)

The "lesson" here is (considering that this was the problem), faced with an error, test your code line by line and observe the result of each step, often you will find the source of error so.