MLE in R language

Question

I am trying to estimate parameters of Birnbaum-Saunders distribution by using its log-likelihood function and 'nlminb' function with R language. But the result differs from the expected one.

My code simply generates sample of BS distribution with parameters alpha=1 and beta = 0.5, then I try to find these values by maximizing the log-likelihood function, but estimated values are "$par 1 0.7959565 0.0001000", which seems quite different from the original one. As I am expecting, these values should converge to the original parameters with increasing of 'n', but in my case it does not work. Here is my code:

n <- 10000
a <- 1
b <- 0.5
sample <- NULL
for(i in 1:n)
{
  z <- rnorm(1, 0, 1)
  t <- b/4 * (a*z + sqrt(4 + (a*z)^2))^2
  sample[i] <- t
}
hist(sample, freq=FALSE)

LLbs <- function(params)
{
  y <- -n*log(params[1]) -n*log(params[2]) +
  sum( log((params[2]/sample[i])^(1/2)+(params[2]/sample[i])^(3/2)))-
  sum( (sample[i]/params[2])+(params[2]/sample[i])-2)/(2*params[1]^2)
  return(-y)
}
nlminb(c(1, 1), LLbs, lower = 0.0001, upper = Inf)

To generate sample I use value of variable T which depends on variable Z which is following standart normal distributions. This link and log-likelihood I found in this research paper (pages 7 and 14).

Answer 1

Remove the [i] from sample[i] and it will work.

This also works, discussed in the paper.

s=mean(sample)
r=1/mean(1/sample)

b=function(b) {
  K=1/mean(1/(b + sample))
  return((b^2 - b*(2*r+K) + r*(s+K))^2)
}

b_mle=optim(1, b, method="BFGS")$par

sqrt(s/b_mle + b_mle/r -2)