nloptr optimization - "objective in x0 Returns NA" error

Question

I’m having trouble using the function nloptr of R. I need the optimal value to be used in the function. The optimization will be iterative, using another value that is also being optimized at each step. As I can’t return a value, I’m also having problems with the function digamma: non-numeric argument to mathematical function.

Follows excerpt from the code in which the nloptr returns the error Error in is.nloptr(ret) : objective in x0 returns NA:

 library(nloptr)

 funcao_f = function(m,p) (-log(p/(m+p)))-(m/(m+p))
 funcao_g = function(m,p) (-1/(m+p))
 funcao_a = function(y,p)  digamma(y+p)-digamma(p)

 E_p= function(Ee_m,p) -(sum(y)/Ee_m)*(funcao_f(Ee_m,p) - Ee_m*funcao_g(Ee_m,p))-sum(funcao_a(y,p))

 y = c(rep(1:4, c(34,10,3,3)))
 n=length(y)
 E_max=1e-5

 soma_y=sum(y)
 Ee_m=c()
 Ee_p=c()
 Ee_p_v=c()#quero valor numerico da sol da otimização
 E_n0=c()

 m= mean(y)          
 p= 0.1                 
 i=0

 while(er>E_max){  
 i=i+1
 er=abs(m-Ee_m[i])
 E_n0[i]=n*dnbinom(0,m,p)/(1-dnbinom(0,m,p))

 Ee_m[i]= soma_y/(n+E_n0[i])
 Ee_p_v[i]= nloptr(x0=p, eval_f=E_p, lb=c(0), ub=c(max(yt)), opts = list("algorithm"="NLOPT_GN_DIRECT","maxeval"=10000,"xtol_rel"=1.0e-20),p=p)
 Ee_p[i]= Ee_p_v[i]$solution[1]#quero valor numerico da sol da otimização

 m=Ee_m[i]
 p=Ee_p_v[i]
 }#fim loop 

Can someone please help me??

Thanks ae!

Answer 1

It is difficult to understand this function. It follows a code that works, but not if you do what you want.

E_p= function(Ee_m,p,y=y) -(sum(y)/Ee_m)*(funcao_f(Ee_m,p) - Ee_m*funcao_g(Ee_m,p))-sum(funcao_a(y,p))
y = c(rep(1:4, c(34,10,3,3)))
y=as.numeric(y)
n=length(y)
E_max=1e-5

soma_y=sum(y)
Ee_m=c()
Ee_p=c()
Ee_p_v=list()#quero valor numerico da sol da otimização
E_n0=c()

m= mean(y)          
p= 0.1                 
i=0
er=1e99
while(er>E_max){  
  i=i+1

  E_n0[i]=n*dnbinom(0,m,p)/(1-dnbinom(0,m,p))

  Ee_m[i]= soma_y/(n+E_n0[i])
  Ee_p_v[[i]]= nloptr(x0=Ee_m[i], eval_f=E_p, lb=c(0), ub=c(max(y)), opts = list("algorithm"="NLOPT_GN_DIRECT","maxeval"=10000,"xtol_rel"=1.0e-20),p=p,y=y)
  Ee_p[i]= Ee_p_v[[i]]$solution[1]#quero valor numerico da sol da otimização
  er=abs(m-Ee_m[i])
  m=Ee_m[i]
  p=Ee_p[i]
}