So I’m interested in finding an approximate maximum entropy distribution for a positive quantity with a given mean and a given mean absolute difference between two randomly chosen values (it’s related to the Gini coefficient).
The analytical solution seems intractable, so I’m thinking to do a numerical approximate solution. I think it’s equivalent for my purposes to have the mean = 1 and the mean absolute difference equal to a given fraction (between 0 and 1 but practically speaking between 0.2 and 0.6)
My thought therefore was to represent the PDF as something like exp(-F(x)) truncated to something like [0,10] where F is a radial basis function expansion, then use a numerical constrained optimizer to maximize the entropy subject to the constraints of normalization, the given mean = 1, and the given Gini coefficient… Calculating the Gini requires a double integral over [0,10] x [0,10] which will have to be done at each step!
There are many Optimization packages available in Julia, and many integration packages. Any suggestions for a constrained optimizer that is likely to handle this well? Any suggestions for a 2D integration method that would handle the Gini calculation well?
So Cubature.jl is working just fine for the 1 and 2D integration I need. What is going to work as a constrained optimizer?
It’s hard to say without seeing your mathematical formulation but have you looked at InfiniteOpt.jl?
I have not. I will look at that. It might be a better choice.
I did manage to make NLopt work for this problem. I’ll come back and post how tomorrow.
1 Like
The way I did this was to represent the pdf as
\exp\left(\sum_{i=1}^N a_i \sqrt{1 + ((x - c_i)/s)^2}\right)/Z
With fixed centers c_i and unknown basis expansion coefficients a_i. I then wrote functions which normalized this, calculated the entropy, and calculated the gini… and used NLopt as follows:
for g in collect(0.25:0.05:0.45)
opt = Opt(:LN_COBYLA,13)
opt.max_objective = (x,g) -> ent(x)
equality_constraint!(opt, (x,gg) -> gin(x) - g,0.005)
equality_constraint!(opt, (x,g) -> muval(x) - 1.0,0.005)
opt.maxtime=120
opt.ftol_rel = 3e-3
opt.xtol_rel = 1e-2
opt.xtol_abs = 1e-3
(optf,optx,ret) = optimize(opt,startx)
startx = optx ./ 2.0
@show(ret)
pp = PdfRbf(centers,optx,3.0,1.0)
normalize!(pp)
push!(pps,pp)
end
It takes tens of seconds to do the calculation, but it’s no big deal, it’s just a one-time calc and then I do my other calculations on the curves.
