# Maximum likelihood for t-distribution by gradient descent

**URL:** <https://discourse.julialang.org/t/maximum-likelihood-for-t-distribution-by-gradient-descent/97954>\
**Category:** Statistics\
**Tags:** question, package\
**Created:** [April 26, 2023, 3:54pm UTC](https://discourse.julialang.org/t/maximum-likelihood-for-t-distribution-by-gradient-descent/97954 "2023-04-26T15:54:05Z")\
**Posts on this page:** 5\
**Page:** 2

<div class="post-metadata">

**Author:** ![kellertuer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kellertuer/32/220707_2.png) [@kellertuer](https://discourse.julialang.org/u/kellertuer)\
**Post date:** [April 27, 2023, 2:06pm UTC](https://discourse.julialang.org/t/maximum-likelihood-for-t-distribution-by-gradient-descent/97954/21 "2023-04-27T14:06:43Z")

</div>

I started fiddling a bit around (first with the Euclidan gradients) with the given formula from the paper above but got a bit stuck.

1. They suddenly take the gradient with respect to \Sigma^{-1}. What is the better variable to use \Sigma or \Sigma^{-1}? Both are SPD anyways.
2. In their gradient with respect tu \nu they suddenly use a function \psi but they either never define that or I missed it (or it is an obvious function in statistics?) _edit:_ Ah I missed it, its the digamma function (luckily available in SpecialFunctions).

_edit2:_ besides the question which variable to take as an optimisation variable, \Sigma or its inverse, I at least managed to implement the three separated gradients. Will continue the next days fiddling around further 🙂

---

<div class="post-metadata">

**Author:** ![dmbates](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dmbates/32/44_2.png) [@dmbates](https://discourse.julialang.org/u/dmbates)\
**Post date:** [April 27, 2023, 2:11pm UTC](https://discourse.julialang.org/t/maximum-likelihood-for-t-distribution-by-gradient-descent/97954/22 "2023-04-27T14:11:13Z")

</div>

Another possibility is the `Distributions.LKJ` or `Distributions.LKJCholesky` parameterization.

Many years ago we used the Cholesky with the log of the diagonal elements in the `nlme` package for R (Pinheiro and Bates, “Mixed-Effects Models in S and S-PLUS”, Springer) to express a positive definite covariance matrix. But in mixed-effects models it is not uncommon to converge to a singular covariance matrix, or need to evaluate the log-likelihood at a singular covariance matrix during the iterative optimization. This causes problems when using the logarithm of the diagonal elements. In the lme4 package for R and in `MixedModels.jl` we switched back to the Cholesky factor itself optimizing with box constraints using Powell’s BOBYQA algorithm as implemented in `NLopt.jl`

The behavior at the boundary can be an important consideration. Often the expressions of the log-likelihood for a mixed-effects model are written in terms of the precision matrix, which is the inverse of the covariance matrix, and many implementations are therefore written in terms of the precision matrix. But the precision doesn’t exist if the covariance matrix is singular, as can and does happen. It can be shown that the other limiting case, a singular precision matrix, can’t be the estimate for these models.

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [April 27, 2023, 3:25pm UTC](https://discourse.julialang.org/t/maximum-likelihood-for-t-distribution-by-gradient-descent/97954/23 "2023-04-27T15:25:42Z")

</div>

I found it easier to read the published version of this paper: [https://www.tandfonline.com/doi/full/10.1080/03610926.2018.1445861](https://www.tandfonline.com/doi/full/10.1080/03610926.2018.1445861)

---

<div class="post-metadata">

**Author:** ![kellertuer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kellertuer/32/220707_2.png) [@kellertuer](https://discourse.julialang.org/u/kellertuer)\
**Post date:** [April 28, 2023, 6:47am UTC](https://discourse.julialang.org/t/maximum-likelihood-for-t-distribution-by-gradient-descent/97954/24 "2023-04-28T06:47:34Z")

</div>

They still have the gradient with the inverse of \Sigma _and_ the typo in Equation (5).

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [April 28, 2023, 5:37pm UTC](https://discourse.julialang.org/t/maximum-likelihood-for-t-distribution-by-gradient-descent/97954/25 "2023-04-28T17:37:38Z")

</div>

Metida.jl work with optimization maximum log-likelihood estimate for MVNormal distributions, see [details](https://pharmcat.github.io/Metida.jl/stable/details/).

You can use gradient descend, Newton or any other method (NelderMead). First of all you should think about constructing Σ matrix in terms of variation and coefficients of variations. Second thoughts - how to make smooth optimization borders: you can use `exp-log` transformation for variation and sigmoid (or arctg) function for `rho` (coefficient of covariation, that can be only in range from -1 to 1).

Inverse of Σ is more performance issue, you can make your own approach if your Σ matrix is diagonal or block-diagonal, also [sweep operator](https://github.com/joshday/SweepOperator.jl) can bu useful.

I made analytical gradient function for log-likelihood, but it was not so faster that gradient with ForardDiff. Original paper can be found [here](https://www.tandfonline.com/doi/abs/10.1080/01621459.1988.10478693).

 ![image](https://global.discourse-cdn.com/julialang/original/3X/8/6/863793b44cbd846480b0a8c20206525937be46f7.png)

Upd:

So, your task about **t-distribution** , and all above can’t be applicable forward, but maybe you can find some insight in this materials.

[Previous page](https://discourse.julialang.org/t/maximum-likelihood-for-t-distribution-by-gradient-descent/97954.md?page=1)