Error in output of the optimum function when creating a regression structure

1

When running the simulation below I get the error at the end of the code, can anyone help me 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

    
asked by anonymous 21.05.2016 / 19:48

1 answer

1

It seems that Eduardo's suggestion in the comments has to do with his problem, yes. In that answer for that question, the tip was given to check the result of your function with the initial values. In your case:

> Loglik(start, t)
[1] NA

This indicates that something is wrong with your function. Investigating it line by line, we can see a problem here:

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

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

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

Notice that the result is now not of NA s:

> alpha
[1] 2 3

And the optim function works 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 you with the evaluation of this result, because you know better what to expect from so much mathematics. :)

The "lesson" here is (considering that this was the problem), facing an error, testing your code line by line and observing the result of each step, you will often find the source error as well.     

22.05.2016 / 00:58