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