# How to fit a Dirichlet Distribution using MLE in Distributions.jl?

**URL:** <https://discourse.julialang.org/t/how-to-fit-a-dirichlet-distribution-using-mle-in-distributions-jl/49645>\
**Category:** Statistics\
**Tags:** statistics, turing, distributions\
**Created:** [November 5, 2020, 9:30pm UTC](https://discourse.julialang.org/t/how-to-fit-a-dirichlet-distribution-using-mle-in-distributions-jl/49645 "2020-11-05T21:30:54Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![stochastic\_guy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stochastic_guy/32/9877_2.png) [@stochastic\_guy](https://discourse.julialang.org/u/stochastic_guy)\
**Post date:** [November 5, 2020, 9:30pm UTC](https://discourse.julialang.org/t/how-to-fit-a-dirichlet-distribution-using-mle-in-distributions-jl/49645/1 "2020-11-05T21:30:54Z")

</div>

I would like to fit Dirichlet Distribution to `all_patterns` using MLE in the Distributions package  
where alpha to Dirichlet is `fill(1,k)` and `all_patterns` is a combination array

```julia
all_patterns = (p=(a,c)->c<2 ? a : [[x;y] for x=a for y=p(a,c-1)])(0:1,K)
all_patterns = convert(Array,VectorOfArray(all_patterns)')

```

and this is what I tried [but not working]

`fit_mle(Dirichlet,fill(1,K),all_patterns)`

using the modified interface specified here: [Distribution Fitting · Distributions.jl](https://juliastats.org/Distributions.jl/stable/fit/)

---

<div class="post-metadata">

**Author:** ![mthelm85](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mthelm85/32/224164_2.png) [@mthelm85](https://discourse.julialang.org/u/mthelm85)\
**Post date:** [November 5, 2020, 11:40pm UTC](https://discourse.julialang.org/t/how-to-fit-a-dirichlet-distribution-using-mle-in-distributions-jl/49645/2 "2020-11-05T23:40:51Z")

</div>

Do you need to use floating-point numbers? The code below works:

```julia
using Distributions

K = 3
all_patterns = (p=(a,c)->c<2 ? a : [[x;y] for x=a for y=p(a,c-1)])(1:2,K)
all_patterns = hcat(all_patterns...)

julia> fit(Dirichlet, Float64.(all_patterns))
Dirichlet{Float64}(alpha=[1.677728917182107e-5, 1.677728917182107e-5, 1.677728917182107e-5])

```

---

<div class="post-metadata">

**Author:** ![stochastic\_guy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stochastic_guy/32/9877_2.png) [@stochastic\_guy](https://discourse.julialang.org/u/stochastic_guy)\
**Post date:** [November 6, 2020, 7:05am UTC](https://discourse.julialang.org/t/how-to-fit-a-dirichlet-distribution-using-mle-in-distributions-jl/49645/3 "2020-11-06T07:05:46Z")

</div>

1. My `all_patterns` is a matrix of binary values (0 & 1) whereas in your code it’s not the case  
1.1 While running you code with binary matrix, I am getting the error  
`ArgumentError: Dirichlet: alpha must be a positive vector.`

2. Dirichlet takes a vector (alphas) as an argument. In your code I think `all_patterns` are being passed as the alphas. I am not sure though.

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [November 6, 2020, 2:24pm UTC](https://discourse.julialang.org/t/how-to-fit-a-dirichlet-distribution-using-mle-in-distributions-jl/49645/4 "2020-11-06T14:24:13Z")

</div>

Did you see my answer on Slack? It would be helpful if you could turn your example into an MWE that people can run, rather than trying to guess what `all_patterns` is (as `K` and `VectorOfArrays` are undefined).

I also don’t understand why in your OP you seem to want to pass a vector of alphas, given that that’s the only parameter of the Dirichlet what are you looking to estimate?

---

<div class="post-metadata">

**Author:** ![stochastic\_guy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stochastic_guy/32/9877_2.png) [@stochastic\_guy](https://discourse.julialang.org/u/stochastic_guy)\
**Post date:** [November 6, 2020, 4:40pm UTC](https://discourse.julialang.org/t/how-to-fit-a-dirichlet-distribution-using-mle-in-distributions-jl/49645/5 "2020-11-06T16:40:08Z")

</div>

Sorry I missed your response in slack. I just saw your reply.

- VectorOfArrays come from here:

> [@Very best way to concatenate an array of arrays](https://discourse.julialang.org/t/very-best-way-to-concatenate-an-array-of-arrays/8672/2):
>
> using RecursiveArrayTools VA = VectorOfArray(A) arr = convert(Array,VA) VectorOfArray is a lazy representation of the array of arrays as a higher order tensor (/matrix if A is a vector of vectors). It uses the indexing fallback to do the conversion. Honestly, I tried writing my own loop for the conversion but the fallback was faster (Tim Holy Cartesian magic?) so this is what I’ve found to be the fastest, and I don’t know how to recreate it without using the lazy matrix indexing.

- Alphas are the priors to the Dirichlet Distribution.
- K should take +ve integer value (User can decide what value to take)

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [November 6, 2020, 7:14pm UTC](https://discourse.julialang.org/t/how-to-fit-a-dirichlet-distribution-using-mle-in-distributions-jl/49645/6 "2020-11-06T19:14:41Z")

</div>

Thanks, that clarifies things a little bit. I still don’t fully understand what you’re trying to do as I haven’t thought through the dimensions of your problem, but it now seems to me that when you say “prior” you mean an initial value for alpha (rather than a prior in the Bayesian sense). In that case you probably need the `init` keyword:

```julia
using Distributions, RecursiveArrayTools

K = 3

all_patterns = (p=(a,c)->c<2 ? a : [[x;y] for x=a for y=p(a,c-1)])(0.0:1.0,K)
all_patterns = convert(Array,VectorOfArray(all_patterns)')

fit_mle(Dirichlet, all_patterns, init = [1. for _ ∈ 1:8])

```

Note that here I used an init vector of length 8, as your `all_patterns` matrix ends up being 8x3 - that’s what I meant by I haven’t thought through the dimensions of your problem, it seems from your initial post that you might want 3x8 here?

---

<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:** [November 7, 2020, 9:53pm UTC](https://discourse.julialang.org/t/how-to-fit-a-dirichlet-distribution-using-mle-in-distributions-jl/49645/7 "2020-11-07T21:53:35Z")

</div>

You could just write a call to Optim.jl to optimize the logpdf.