Как я могу оценить структурный параметр модели поиска работы с помощью otpim

Question

Мне сложно оценить параметр моей модели. Это структурная модель поиска работы и соответствия. Я пишу, используя Optim и BFGS метод. Но многие предупреждения появляются как Предупреждающие сообщения:

1: In pnorm(log(wstar), mean = mu, sd = sigma, log.p = FALSE) : NaNs produced
2: In pnorm(log(wstar), mean = muH, sd = sigmaH, log.p = FALSE) : NaNs produced
3: In pnorm(log(wstar + d), mean = muH, sd = sigmaH, log.p = FALSE) :
  NaNs produced 

Мой код

    ll<-function(par){
  lamda  <-par[1]
  lamdaH <-par[2]
  delta  <-par[3]
  deltaH <-par[4]
  mu     <-par[5]
  muH    <-par[6]
  sigma  <-par[7]
  sigmaH <-par[8]
  pi     <-par[9]
  d      <-par[10]
  wstar<- min(Data$wh[Data$wh!=0],na.rm=TRUE)
  NeA<-table(Data$handicap[Data$status==1])[[1]]
  NeD<-table(Data$handicap[Data$status==1])[[2]]
  NuA<-table(Data$handicap[Data$status==2])[[1]]
  NuD<-table(Data$handicap[Data$status==2])[[2]]
  h<-lamda*(1- pnorm(log(wstar), mean = mu, sd = sigma,log.p=FALSE))
  hh<-(lamdaH*((1-pi)*(1- pnorm(log(wstar), mean = muH, sd = sigmaH,log.p=FALSE))
             + pi*(1- pnorm(log(wstar+d), mean = muH, sd = sigmaH,log.p=FALSE))))
  wn <-(1/alpha)*(Data$wh[Data$wh>0]-(1-alpha)*wstar)
  wp <-(1/alpha)*(Data$wh[Data$wh>0]-(1-alpha)*wstar +alpha*d)

  fe <-(1/alpha*(dnorm(log(na.omit(wn[wn>0])),mean=mu,sd=sigma, log=FALSE)
             /(sigma*na.omit(wn[wn>0]))/(1-pnorm(log(wstar),mean = mu,sd=sigma,log.p=FALSE))))
  fn <-((1-pi)/alpha*(dnorm(log(na.omit(wn[wn>0])),mean=muH,sd=sigmaH, log=FALSE)
                   /(sigmaH*na.omit(wn[wn>0]))/(1-pnorm(log(wstar),mean = muH,sd=sigmaH,log.p=FALSE))))
  fp <-(pi/alpha*(dnorm(log(na.omit(wp[wp>0])),mean=muH,sd=sigmaH, log=FALSE)
               /(sigmaH*na.omit(wp[wp>0]))/(1-pnorm(log(wstar+d),mean = muH,sd=sigmaH,log.p=FALSE))))

  llh<-((NuA+NeA)*log(h/(delta+h))-h*sum(na.omit(Data$spellU[Data$handicap==0]))+NuA*log(delta)
  +sum(na.omit(log(fe)[Data$status==1&Data$handicap==0]))
  +(NuD+NeD)*log(hh/(deltaH+hh))-hh*sum(na.omit(Data$spellU[Data$handicap==1]))+NuD*log(deltaH)
  +sum(na.omit(log((fn+fp)[Data$status==1&Data$handicap==1]))))

  if(!is.finite(llh))
  llh <-1e+20 
  return(llh)}
para     <-c(0.20,0.10,0.009,0.005,4.3,3.2,0.52,0.32, 0.52,10)
res     <- optim(par=para,fn=ll,method="BFGS", gr = function(Data) pracma::grad(ll, Data),control = list(factr=1e-10,trace = TRUE,REPORT = 1,maxit= 10000), hessian=TRUE)