I'm trying to find an optimal solution for g40 and g60, such that change=0.01 Below is my code, however i'm getting the error Error in .check_function_for_sanity(F, n) : cannot evaluate function 'F' using 'n' = 2 parameters
I tried to look for a solution and I found that I should define a wrapper function, but it still doesn't work. I'm new for optimization on R, so if there is another mistake, do not hesitate to teach me.
Below is my code (the data file contains two columns: p and L)
#weights
p40=which(abs(data$p-0.4)==min(abs(data$p-0.4))) #position of the 40th percentile
w40=data[p40,2]
w60=1-w40
#Non-cumulative L
L_non=data$L[1]
for(i in 2:nrow(data))
L_non[i]=data$L[i]-data$L[i-1]
#X
x1=sum(data$p*L_non)/sum(L_non)
x2=sum(data$p[1:p40]*L_non[1:p40])/sum(L_non[1:p40])
#Define function for the optimisation
FG=function(g40,g60) {
#Gammas and m
gama=g40*w40+g60*w60
gama40=g40
m=gama40-gama
#theta & delta
theta=-m/(x2-x1)
delta=gama+theta*x1
#L_prime
L_prime=(1+delta-theta*data$p)*L_non
data=cbind(data,L_prime)
#Gini approximated from the new distribution
Lp0.2=data[which(abs(data$p-0.2)==min(abs(data$p-0.2))),3]
Lp0.4=data[which(abs(data$p-0.4)==min(abs(data$p-0.4))),3]
Lp0.6=data[which(abs(data$p-0.6)==min(abs(data$p-0.6))),3]
Lp0.8=data[which(abs(data$p-0.8)==min(abs(data$p-0.8))),3]
Gini_prime=2*(125/288)*(1/5)*(3*(0.2-Lp0.2)+2*(0.4-Lp0.4)+2*(0.6-Lp0.6)+3*(0.8-Lp0.8))
#change
change=(Gini_prime/Gini)-1
}
# Check that the F_objective function works:
wrapper <- function(x) FG(x, 0.1)
#OP problem
library(ROI)
library(ROI.plugin.nloptr)
prob<-OP(objective=F_objective(F = wrapper, n = 2),
constraints=F_constraint(sum(L_prime), dir=c("=",),rhs=c(1)) & L_constraint(change, dir=c("=",),rhs=c(0.01))
)
solve=ROI_solve(prob,solver="nloptr.cobyla",start=c(0,0))
H=solution(solve)
#weights
p40=which(abs(data$p-0.4)==min(abs(data$p-0.4))) #position of the 40th percentile
w40=data[p40,2]
w60=1-w40
#Non-cumulative L
L_non=data$L[1]
for(i in 2:nrow(data))
L_non[i]=data$L[i]-data$L[i-1]
#X
x1=sum(data$p*L_non)/sum(L_non)
x2=sum(data$p[1:p40]*L_non[1:p40])/sum(L_non[1:p40])
################################### Define Optimization Functions ###################################
FG=function(g) { # Objective function
#Gammas and m
gama=g[1]*w40+g[2]*w60
gama40=g[1]
m=gama40-gama
return(m)
}
FP=function(g) { #sum of L_prime
#Gammas and m
gama=g[1]*w40+g[2]*w60
gama40=g[1]
m=gama40-gama
#theta & delta
theta=-m/(x2-x1)
delta=gama+theta*x1
#L_prime
L_prime=(1+delta-theta*data$p)*L_non
return(sum(L_prime))
}
CH=function(g){ #Gini calculator
#Gammas and m
gama=g[1]*w40+g[2]*w60
gama40=g[1]
m=gama40-gama
#theta & delta
theta=-m/(x2-x1)
delta=gama+theta*x1
#L_prime
L_prime=(1+delta-theta*data$p)*L_non
#Cumulative L
LC_prime=NULL
LC_prime[1]=L_prime[1]
for(i in 2:length(L_prime)) LC_prime[i]=LC_prime[i-1]+L_prime[i]
#Gini approximated from the new distribution
Lp0.2=LC_prime[which(abs(data$p-0.2)==min(abs(data$p-0.2)))]
Lp0.4=LC_prime[which(abs(data$p-0.4)==min(abs(data$p-0.4)))]
Lp0.6=LC_prime[which(abs(data$p-0.6)==min(abs(data$p-0.6)))]
Lp0.8=LC_prime[which(abs(data$p-0.8)==min(abs(data$p-0.8)))]
Gini_prime=2*(125/288)*(1/5)*(3*(0.2-Lp0.2)+2*(0.4-Lp0.4)+2*(0.6-Lp0.6)+3*(0.8-Lp0.8))
#change
# change=(Gini_prime/Gini)-1
return(Gini_prime)
}
################################### OP problem ###################################
library(ROI)
library(ROI.plugin.nloptr)
prob<-OP(objective=F_objective(F = FG, n = 2),
constraints=F_constraint(list(F=CH,F=FP), dir=c("<=","<="),rhs=c(Gini*growth,1)),
types=c("C","C"),
bounds = V_bound(ui = seq_len(2), lb = c(-5,-5),ub=c(5,5))
)
solve=ROI_solve(prob,solver="nloptr.cobyla",start=c(0,0))
G=solution(solve)