I created the code below to parameter estimation of the Birnbaum-Saunders distribution, which most want is the implementation of the gradient, which was added to the code! I wonder if it is correct? And if there is a better way to automate the gradient with others any packages, for example the NumDeriv, Deriv or others?rm(list=ls())cat(\014)library(VGAM)library(Deriv)n=100N=30m=matrix(ncol=2,nrow=N)f <- expression((log((1/(2*alpha*beta*sqrt(2*pi)))*(sqrt(beta/t)+(beta/t)^(3/2))*exp(-(1/(2*alpha^2))*((t/beta)+(beta/t)-2)))))gd1 <- Deriv(f,alpha)gd2 <- Deriv(f,beta)for (i in 1:N){  alpha<-2  beta <-1x <- rnorm(n,mean = 0,sd=sqrt((alpha^2)/4))t <- beta*(1+2*x^2+2*x*sqrt(1+x^2)) #t~BS(alpha,beta)#t <- rbisa(n, alpha, beta)#sum(1*(t<0))#theta <- c(1,2)rm(alpha,beta)log.ver <- function(theta){  alpha <- theta[1]  beta  <- theta[2]  l <- sum(log((1/(2*alpha*beta*sqrt(2*pi)))*(sqrt(beta/t)+(beta/t)^(3/2))*exp(-(1/(2*alpha^2))*((t/beta)+(beta/t)-2))))  return(-l)} grad<-function(theta){  alpha <- theta[1]  beta  <- theta[2]  gd <- cbind(eval(gd1),eval(gd2))  return(-colSums(gd, na.rm = FALSE))  #return(c(-eval(gd1),-eval(gd2)))  }alpha_0 <- 3beta_0 <- 2start <- c(alpha_0,beta_0)opt <- optim(start,log.ver,method=BFGS,gr=grad,hessian = F)m[i,]=opt$par}alpha<-2beta <-1#Calculating the average of each column of the array of parameters mmest=colMeans(m)#calculating the standard deviation of each column of the array of parameters mdest=apply(m,2,sd)#root mean square error in the calculation of each column of the array of parameters m in relation to the true value of the parametereqm=function(x,opt){   N=length(x)  sqrt(sum(((x-opt)^2))/N)}#Estimated mean squared error of each parameter eqmest=c(eqm(x=m[,1],opt=alpha),         eqm(x=m[,2],opt=beta))# Table with the true values of the parameters and the average# Standard deviation and mean square error of the estimated parametersmtruevalue <- c(alpha,beta)tab=data.frame(truevalue,mean=mest,sd=dest,eqm=eqmest)tabMaybe there is a better/faster more elegant way to do it, if yes please share! Thanks!!
gradient vector for use in the optim function
r;simulation