# ModelingToolkit & Distributions for MLE estimation

**URL:** <https://discourse.julialang.org/t/modelingtoolkit-distributions-for-mle-estimation/97613>\
**Category:** Optimization (Mathematical)\
**Tags:** optimization, distributions, mle, modelingtoolkit, symbolics\
**Created:** [April 18, 2023, 12:33pm UTC](https://discourse.julialang.org/t/modelingtoolkit-distributions-for-mle-estimation/97613 "2023-04-18T12:33:06Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![dmetivie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dmetivie/32/6926_2.png) [@dmetivie](https://discourse.julialang.org/u/dmetivie)\
**Post date:** [April 18, 2023, 12:33pm UTC](https://discourse.julialang.org/t/modelingtoolkit-distributions-for-mle-estimation/97613/1 "2023-04-18T12:33:07Z")

</div>

I am trying to use [ModelingToolKit.jl](https://github.com/SciML/ModelingToolkit.jl) to write a maximum loglikelihood estimator for generic distributions eg Univariate, Mixture, Multivariate etc.

```julia
using ModelingToolkit, Optimization, OptimizationOptimJL
using Distributions

f(θ) = MixtureModel([Exponential(θ[1]), Exponential(θ[2])], [θ[3], 1-θ[3]])
ℓ(θ, x) = -loglikelihood(f(θ), x) # = sum(logpdf(f(θ), y) for y in x)

# Data
N = 100
θtrue = [10, 1, 0.2]
x = rand(f(θtrue), N)

# Observation as parameters (not sure if this is needed?)
@parameters y[1:N]

# variable
@variables θ[1:3]

mle = ℓ(θ, y)

# # Not tested yet, the code above need to work first
# @named sys = OptimizationSystem(mle, θ, y)
# θ0 = [θ[1] => 20, θ[2] => 4, θ[3] => 0.5]
# p = [y[i] => x[i] for i in eachindex(x)]

# prob = OptimizationProblem(sys, θ0, grad = true, hess = true)

```

However, this does not work:

```julia
ERROR: TypeError: non-boolean (Num) used in boolean context
....
[4] Exponential(θ::Num)
@ Distributions C:\Users\dmetivie\.julia\packages\Distributions\Spcmv\src\univariate\continuous\exponential.jl:29

```

I tried `@register_symbolic Exponential(θ)` but now the error is

```julia
ERROR: MethodError: no method matching MixtureModel(::Vector{Num}, ::Vector{Num})
Closest candidates are:
  MixtureModel(::Type{C}, ::AbstractArray) where C<:Distribution at C:\Users\metivier\.julia\packages\Distributions\Spcmv\src\mixtures\mixturemodel.jl:127
  MixtureModel(::Type{C}, ::AbstractArray, ::Vector{T}) where {C<:Distribution, T<:Real} at C:\Users\metivier\.julia\packages\Distributions\Spcmv\src\mixtures\mixturemodel.jl:144
  MixtureModel(::Vector{C}, ::VT) where {C<:Distribution, VT<:(AbstractVector{<:Real})} at C:\Users\metivier\.julia\packages\Distributions\Spcmv\src\mixtures\mixturemodel.jl:138
Stacktrace:
 [1] f(θ::Symbolics.Arr{Num, 1})
   @ Main c:\Users\metivier\Dropbox\PC (2)\Documents\Simulations\julia_mwe\ModelingToolkit_test_for_TrigMixture\mwe_modelingtoolkit_discourse.jl:4
 [2] ℓ(θ::Symbolics.Arr{Num, 1}, x::Symbolics.Arr{Num, 1})
   @ Main c:\Users\metivier\Dropbox\PC (2)\Documents\Simulations\julia_mwe\ModelingToolkit_test_for_TrigMixture\mwe_modelingtoolkit_discourse.jl:5
 [3] top-level scope
   @ c:\Users\metivier\Dropbox\PC (2)\Documents\Simulations\julia_mwe\ModelingToolkit_test_for_TrigMixture\mwe_modelingtoolkit_discourse.jl:25

```

It seems that `Distributions.jl` and `ModelingToolKit.jl` are not compatible in general, since there are type check, `logsumexp` with if statements, etc do not like `Num`.  
I saw that [DistributionsAD.jl](https://github.com/TuringLang/DistributionsAD.jl) was a good candidate to answer the issue. However, I was not able to understand how it works (could not find documentation).

Anyone with an idea how to implement that?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 18, 2023, 2:30pm UTC](https://discourse.julialang.org/t/modelingtoolkit-distributions-for-mle-estimation/97613/2 "2023-04-18T14:30:15Z")

</div>

Why do you need to use ModelingToolkit here, as opposed to just calling Optimization.jl directly?

---

<div class="post-metadata">

**Author:** ![dmetivie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dmetivie/32/6926_2.png) [@dmetivie](https://discourse.julialang.org/u/dmetivie)\
**Post date:** [April 18, 2023, 3:15pm UTC](https://discourse.julialang.org/t/modelingtoolkit-distributions-for-mle-estimation/97613/3 "2023-04-18T15:15:57Z")

</div>

For reference this works

```julia
using Optimization, OptimizationOptimJL
using Distributions

f(θ) = MixtureModel([Exponential(θ[1]), Exponential(θ[2])], [θ[3], 1-θ[3]])
ℓ(θ, x) = -loglikelihood(f(θ), x) # = sum(logpdf(f(θ), y) for y in x)

# Data
N = 100
θtrue = [10, 1, 0.2]
x = rand(f(θtrue), N)

θ0 = [20, 4, 0.5]
_p = x

OptFunc = OptimizationFunction(ℓ, Optimization.AutoForwardDiff())
l1 = ℓ(θ0, _p)
# prob = OptimizationProblem(OptFunc, θ0, _p) # it goes out of legal distributions bound!
prob = OptimizationProblem(OptFunc, θ0, _p, lb = [0.01, 0.01, 0.01], ub = [100, 100, 1])

# Start with some derivative-free optimizers
sol = solve(prob, SAMIN())

```

Does it mean that the problem is AutomaticDifferentiation compliant ?

Anyway, I wanted to use ModelingToolKit to get some of what they advertise on the [doc](https://docs.sciml.ai/ModelingToolkit/stable/tutorials/optimization/#Nested-Systems) "It provides many other features like auto-parallelism and sparsification too. Plus, you can hierarchically nest systems to generate huge optimization problems. ".

Would it make a difference here? I don’t know (guess not).  
However, for the true problem I have in mind, the simplications of ModelingToolKit should be very useful.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 19, 2023, 7:31am UTC](https://discourse.julialang.org/t/modelingtoolkit-distributions-for-mle-estimation/97613/4 "2023-04-19T07:31:42Z")

</div>

Its simplifications can help quite a bit on constrained optimization problems and with MOI, but I’m not sure this case would get any benefits (at least with the current integration).

---

<div class="post-metadata">

**Author:** ![dmetivie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dmetivie/32/6926_2.png) [@dmetivie](https://discourse.julialang.org/u/dmetivie)\
**Post date:** [April 19, 2023, 4:17pm UTC](https://discourse.julialang.org/t/modelingtoolkit-distributions-for-mle-estimation/97613/5 "2023-04-19T16:17:04Z")

</div>

Sure, here it is a MWE to try to understand what is not working.  
In my real application, I have periodic terms in the loglikelihood, analytically one can factor all terms with same value to optimize a shorter function.  
Something like

```julia
θ[1]*y[1] + θ[2]*y[2] + ... + θ[1]*y[100] + θ[2]*y[101] + ... = θ[1]*sum(y[i] ...) + θ[2]*sum(y[i] ...)

```

That is where I believed ModelingToolKit could shine.

However, that little experiment triggered questions:

- it seems that all `logpdf` I use are AD compatible with Optimisation (even `logpdf` of mixture). I believed it was not the case (that is why `DistributionsAD` existed). Am I missing something?
- I read somewhere that if it was AD compatible, it was basically compatible with symbolic packages like ModelingToolKit. Here it does not seem to be true.