# Numerical Maximum Entropy... suggested method / Optimization & Integration packages?

**URL:** <https://discourse.julialang.org/t/numerical-maximum-entropy-suggested-method-optimization-integration-packages/85367>\
**Category:** General Usage\
**Tags:** optimization\
**Created:** [August 5, 2022, 7:57pm UTC](https://discourse.julialang.org/t/numerical-maximum-entropy-suggested-method-optimization-integration-packages/85367 "2022-08-05T19:57:10Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [August 5, 2022, 7:57pm UTC](https://discourse.julialang.org/t/numerical-maximum-entropy-suggested-method-optimization-integration-packages/85367/1 "2022-08-05T19:57:10Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [August 5, 2022, 11:37pm UTC](https://discourse.julialang.org/t/numerical-maximum-entropy-suggested-method-optimization-integration-packages/85367/2 "2022-08-05T23:37:49Z")

</div>

So Cubature.jl is working just fine for the 1 and 2D integration I need. What is going to work as a constrained optimizer?

---

<div class="post-metadata">

**Author:** ![jd-foster](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jd-foster/32/35824_2.png) [@jd-foster](https://discourse.julialang.org/u/jd-foster)\
**Post date:** [August 6, 2022, 1:13am UTC](https://discourse.julialang.org/t/numerical-maximum-entropy-suggested-method-optimization-integration-packages/85367/3 "2022-08-06T01:13:57Z")

</div>

It’s hard to say without seeing your mathematical formulation but have you looked at InfiniteOpt.jl?

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [August 6, 2022, 3:45am UTC](https://discourse.julialang.org/t/numerical-maximum-entropy-suggested-method-optimization-integration-packages/85367/4 "2022-08-06T03:45:33Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [August 6, 2022, 5:26pm UTC](https://discourse.julialang.org/t/numerical-maximum-entropy-suggested-method-optimization-integration-packages/85367/5 "2022-08-06T17:26:22Z")

</div>

![image](https://global.discourse-cdn.com/julialang/original/3X/3/c/3cc05883ac9910c9c0e8a8923474b7b4ef0f55cb.png)

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:

```julia
    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.