I wrote a code for a simulation study to obtain average parameter estimates, their mean biases and mean square errors as sample size increases. The r code for the simulation study is presented below.

rm(list=ls())
library(stats4)
library(bbmle)
library(stats)
library(numDeriv)
library(Matrix)

##define PrWG pdf
PrWG_pdf<-function(theta,lambda,gamma,p,x){
  (1-p)*(theta^2*x^(theta-1)+lambda*gamma*x^(gamma-1)*(1+x^theta))*(1+x^theta)^(-theta-1)*exp(-lambda*x^gamma)/(1-(p*(1+x^theta)^(-theta)*exp(-lambda*x^gamma)))^2
}

##define PrWG Quantiles
PrWG_quantile<-function(theta,lambda,gamma,p,u){
  f<-function(x){
    theta*log(1+x^theta)+lambda*x^gamma+log((1-u)/(1-p*u))
  }
  x=uniroot(f,c(0,100),tol=0.0001)$root
  return(x)
}

##Maximum likelihood estimation
PrWG_LL<-function(par){(1-par[4])*(par[1]^2*x^(par[1]-1)+par[2]*par[3]*x^(par[3]-1)*(1+x^par[1]))*(1+x^par[1])^(-par[1]-1)*exp(-par[2]*x^par[3])/(1-(par[4]*(1+x^par[1])^(-par[1])*exp(-par[2]*x^par[3])))^2
}

##Monte-Carlo simulation study
theta=1.0
lambda=0.5
gamma=0.6
p=0.5

nl=c(25,50,100,200,400,800)
for (j in 1:length(nl)){
  n=nl[j]
  N=1000
  mle_theta<-c(rep(0,N))
  mle_lambda<-c(rep(0,N))
  mle_gamma<-c(rep(0,N))
  mle_p<-c(rep(0,N))
  LC_theta<-c(rep(0,N))
  UC_theta<-c(rep(0,N))
  LC_lambda<-c(rep(0,N))
  UC_lambda<-c(rep(0,N))
  LC_gamma<-c(rep(0,N))
  UC_gamma<-c(rep(0,N))
  LC_p<-c(rep(0,N))
  UC_p<-c(rep(0,N))
  
  count_theta=0
  count_lambda=0
  count_gamma=0
  count_p=0
  
  temp=1
  HH1<-matrix(c(rep(2,16)),nrow=4,ncol=4)
  HH2<-matrix(c(rep(2,16)),nrow=4,ncol=4)
  for (i in 1:N)
  {
    print(i)
    flush.console()
    repeat{
      x<-c(rep(0,n))
      
    ##Generate a random variable from uniform distribution
      u<-0
      u<-runif(n,min=0,max=1)
      
      for (m in 1:n){
        x[m]<-PrWG_quantile(theta,lambda,gamma,p,u[m])
      }
      
      mle.result<-nlminb(c(theta,lambda,gamma,p),PrWG_LL,lower=c(0,0,0,0),upper=c(Inf,Inf,Inf,1))
      
      temp=mle.result$convergence
      if(temp==0){
        temp_theta<-mle.result$par[1]
        temp_lambda<-mle.result$par[2]
        temp_gamma<-mle.result$par[3]
        temp_p<-mle.result$par[4]
        HH1<-hessian(PrWG_LL,c(temp_theta,temp_lambda,temp_gamma,temp_p))
        if(sum(is.nan(HH1))==0 & (diag(HH1)[1]>0) & (diag(HH1)[2]>0) & (diag(HH1)[3]>0) & (diag(HH1)[4]>0)){
          HH2<-solve(HH1)
          #print(det(HH1))
        }
        else{
          temp=1}
      }
      if((temp==0) & (diag(HH2)[1]>0) & (diag(HH2)[2]>0) & (diag(HH2)[3]>0) & (diag(HH2)[4]>0) & (sum(is.nan(HH2))==0)){
        break
      }
      else{
        temp=1}
    }
    temp=1
    mle_theta[i]<-mle.result&par[1]
    mle_lambda[i]<-mle.result&par[2]
    mle_gamma[i]<-mle.result&par[3]
    mle_p[i]<-mle.result&par[4]
    
    HH<-hessian(PrWG_LL,c(mle_theta[i],mle_lambda[i],mle_gamma[i],mle_p[i]))
    H<-solve(H)
    LC_theta[i]<-mle_theta[i]-qnorm(0.975)*sqrt(diag(H)[1])
    UC_theta[i]<-mle_theta[i]+qnorm(0.975)*sqrt(diag(H)[1])
    if ((LC_theta[i]<=theta) & (theta<=UC_theta[i])){
      count_theta=count_theta+1
    }
    LC_lambda[i]<-mle_lambda[i]-qnorm(0.975)*sqrt(diag(H)[2])
    UC_lambda[i]<-mle_lambda[i]+qnorm(0.975)*sqrt(diag(H)[2])
    if ((LC_lambda[i]<=lambda) & (lambda<=UC_lambda[i])){
      count_lambda=count_lambda+1
    }
    LC_gamma[i]<-mle_gamma[i]-qnorm(0.975)*sqrt(diag(H)[3])
    UC_gamma[i]<-mle_gamma[i]+qnorm(0.975)*sqrt(diag(H)[3])
    if ((LC_gamma[i]<=gamma) & (gamma<=UC_gamma[i])){
      count_gamma=count_gamma+1
    }
    LC_p[i]<-mle_p[i]-qnorm(0.975)*sqrt(diag(H)[4])
    UC_p[i]<-mle_p[i]+qnorm(0.975)*sqrt(diag(H)[4])
    if ((LC_p[i]<=p) & (p<=UC_p[i])){
      count_p=count_p+1
    }
    
    ##Calculate Average Bias
    ABias_theta<-sum(mle_theta-theta)/N
    ABias_lambda<-sum(mle_lambda-lambda)/N
    ABias_gamma<-sum(mle_gamma-gamma)/N
    ABias_p<-sum(mle_p-p)/N
    print(cbind(ABias_theta,ABias_lambda,ABias_gamma,ABias_p))
    
    ##Calculate RMSE
    RMSE_theta<-sqrt(sum((theta-mle_theta)^2)/N)
    RMSE_lambda<-sqrt(sum((lambda-mle_lambda)^2)/N)
    RMSE_gamma<-sqrt(sum((gamma-mle_gamma)^2)/N)
    RMSE_p<-sqrt(sum((p-mle_p)^2)/N)
    print(cbind(RMSE_theta,RMSE_lambda,RMSE_gamma,RMSE_p))
    
    ##Converge probability
    CP_theta<-count_theta/N
    CP_lambda<-count_lambda/N
    CP_gamma<-count_gamma/N
    CP_p<-count_p/N
    print(cbind(CP_theta,CP_lambda,CP_gamma,CP_p))
    
    ##Average Width
    AW_theta<-sum(abs(UC_theta-LC_theta))/N
    AW_lambda<-sum(abs(UC_lambda-LC_lambda))/N
    AW_gamma<-sum(abs(UC_gamma-LC_gamma))/N
    AW_p<-sum(abs(UC_p-LC_p))/N
    print(cbind(AW_theta,AW_lambda,AW_gamma,AW_p))
  }
}

After running the code in r, the response I got is

Error in hessian.default(PrWG_LL, c(temp_theta, temp_lambda, temp_gamma,: Richardson method for hessian assumes a scalar valued function.

How do I stop the error and get the solution?

Thank you.