2b:null 2c:["$","title",null,{"children":"Distributions on the unit simplex | Mattias Villani | Observable"}] 31:T573,The two main distributions on the unit simplex (${tex`\sum_{k=1}^K x_k = 1 \text{ and } x_k \geq 0`}) are the [Dirichlet distribution](https://observablehq.com/@mattiasvillani/dirichlet-distribution) and the [multivariate LogitNormal](https://observablehq.com/@mattiasvillani/multivariate-logitnormal-distribution) distribution. The Dirichlet distribution generalizes the [Beta distribution](https://observablehq.com/@mattiasvillani/beta-distribution) to the multivariate setting. It plays an important role in Bayesian analysis as the conjugate prior for the multinomial model. It has many tractable properties (moments, marginal distributions etc), but can be a restrictive in that the user has essentially no control over the correlation. The LogitNormal distribution generalizes the univariate [LogitNormal distribution](https://observablehq.com/@mattiasvillani/logit-normal-distribution) to the multivariate setting. The LogitNormal is more flexible than the Dirichlet distribution (thought they are never exactly the same), but many of its properties are not available in closed form. This notebook compares the Dirichlet distribution side-by-side. It also allows you to plot the closest multivariate LogitNormal distribution approximation to the Dirichlet distribution (in the Kullback-Leibler sense). The Ternary plot is from **Herb Susmann** at https://observablehq.com/@herbps1032:T4a2,

Multivariate LogitNormal distribution
The distribution is obtained from a drawing a ${tex`K-1`}-dimensional multivariate normal vector
${tex.block`\boldsymbol{y} =(y_1,\ldots,y_{K-1})^\top \sim N(\boldsymbol{\mu},\boldsymbol{\Sigma})`} and applying a softmax transformation ${tex.block`x_k = \frac{e^{y_k}}{\sum_{k=1}^{K}e^{y_k}} \text{ for }k=1,\ldots,K, \text{ where }y_K = 0`}
The density is given by ${tex.block` \mathrm{LogitNormal}(\boldsymbol{x} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma}) = \vert 2\pi \boldsymbol{\Sigma} \vert^{-1/2} \Big(\prod_{k=1}^K x_k\Big)^{-1}\exp\Big( -\frac{1}{2}(\tilde \boldsymbol{x}-\boldsymbol{\mu})^\top\boldsymbol{\Sigma}^{-1}(\tilde \boldsymbol{x} - \boldsymbol{\mu}) \Big) `} where ${tex`\tilde \boldsymbol{x} = (\log(x_1/x_K),\ldots,\log(x_{K-1}/x_K))^\top`}.

We will here consider the case with ${tex`K=3`} where ${tex.block` \begin{aligned} \boldsymbol{\mu} &= (\mu_1, \mu_2)^\top \\ \boldsymbol{\Sigma} &= \begin{pmatrix} \sigma_1^2 & \rho \sigma_1 \sigma_2\\ \rho \sigma_1 \sigma_2 & \sigma_2^2 \end{pmatrix} \end{aligned} `}33:[] 2f:{"id":"4c0ca43f8c3fea86","slug":"distributions-on-the-unit-simplex","title":"Distributions on the unit simplex","thumbnail":"731309b86f99ecd45949aaa2689a6a5e04b26809b0d84f358fe1e5738f72f08c","access_level":"public","create_time":"$D2023-12-31T10:26:19.958Z","update_time":"$D2026-08-22T18:38:49.490Z","view_time":null,"first_public_time":"$D2023-12-31T11:43:49.136Z","fork":{"id":"14554c2f6979d078","slug":"multivariate-logitnormal-distribution","title":"Multivariate LogitNormal Distribution","owner":{"id":"7592a7ab121f74a0","login":"mattiasvillani","name":"Mattias Villani"}},"stars":1,"starred":false,"owner":{"type":"team","id":"7592a7ab121f74a0","login":"mattiasvillani","name":"Mattias Villani","avatar_url":"https://avatars-next.observableusercontent.com/a440aaf39a3abefb86cbcc758fd7a5ab42555ae68d6825b52ec006f385a94d1d","tier":"pro_2024"},"version":1187,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"56499be8cf43f0ef","login":"mattiasvillani","name":"Mattias Villani","avatar_url":"https://avatars-next.observableusercontent.com/a440aaf39a3abefb86cbcc758fd7a5ab42555ae68d6825b52ec006f385a94d1d","tier":null}],"members":[{"type":"individual","id":"56499be8cf43f0ef","login":"mattiasvillani","name":"Mattias Villani","avatar_url":"https://avatars-next.observableusercontent.com/a440aaf39a3abefb86cbcc758fd7a5ab42555ae68d6825b52ec006f385a94d1d","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"subscribed":null,"forks":0,"collection":{"id":"9fdb7c9ba7be6ba8","slug":"bayesian-learning","title":"Bayesian Learning","owner":{"id":"7592a7ab121f74a0","login":"mattiasvillani","name":"Mattias Villani"}},"type":"notebook","body":{"title":"Distributions on the unit simplex","cells":[{"id":706,"value":"# Distributions on the unit simplex","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":858,"value":"$31","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1041,"value":"

\nDirichlet distribution
\n\nThe distribution is obtained from a drawing a ${tex`K`}-dimensional vector ${tex`\\boldsymbol{y} =(y_1,\\ldots,y_K)^\\top`} with independent elements ${tex`y_k \\sim \\mathrm{Gamma}(\\alpha_k,1)`} and normalizing\n\n${tex.block`x_k = \n\\frac{y_k}{\\sum_{k=1}^K y_k} \\text{ for }k=1,\\ldots,K.`}\n\nThe density is given by\n${tex.block`\n \\mathrm{Dirichlet}(\\boldsymbol{x} \\mid \\boldsymbol{\\alpha}) = \\frac{1}{\\mathrm{B}(\\boldsymbol{\\alpha})} \\prod_{i=1}^k x_i^{\\alpha_i - 1}\n`}\nwhere ${tex`\\mathrm{B}(\\boldsymbol{\\alpha})`} is the multivariate Beta function and ${tex`\\boldsymbol{\\alpha} = \\{ \\alpha_1, \\dots, \\alpha_k \\}`} with ${tex`\\alpha_i > 0`}.","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1056,"value":"$32","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":986,"value":"### Interactive comparison","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1013,"value":"Dirichlet","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":893,"value":"viewof inputs = Inputs.form({\n alpha1: Inputs.range([0.01, 10], {value: 2, label: tex`\\alpha_1`, step: 0.01}),\n alpha2: Inputs.range([0.01, 10], {value: 2, label: tex`\\alpha_2`, step: 0.01}),\n alpha3: Inputs.range([0.01, 10], {value: 2, label: tex`\\alpha_3`, step: 0.01}),\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1071,"value":"Multivariate LogitNormal","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1082,"value":"viewof klcloseparam = Inputs.toggle({label: \"Set parameters to approx Dirichlet\", value: false})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1074,"value":"viewof paramlogitn = { \n if (klcloseparam){\n return Inputs.form([Inputs.text({value: \"LogitNormal parameters set to approximate Dirichlet\", label: \"\", disabled: true})])\n }else{\n return Inputs.form({\n mu1: Inputs.range([-1, 1], {value: 1, label: tex`\\mu_1`, step: 0.1}),\n mu2: Inputs.range([-1, 1], {value: 0.5, label: tex`\\mu_2`, step: 0.1}),\n sigma1: Inputs.range([0.1, 2], {value: 0.5, label: tex`\\sigma_1`, step: 0.1}),\n sigma2: Inputs.range([0.1, 2], {value: 1, label: tex`\\sigma_2`, step: 0.1}),\n rho: Inputs.range([-0.99, 0.99], {value: 0.9, label: tex`\\rho`, step: 0.01})\n })\n }\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":940,"value":"\n \n \n \n
\n
\n
\n \n ${tex`\\mathrm{Dirichlet}\\big(\n \\boldsymbol{\\alpha} = (${inputs.alpha1},${inputs.alpha2},${inputs.alpha3})\\big)`}
\n ${dirichletplt}
\n ${dirichletColorScale}\n
\n
\n \n ${tex`\\mathrm{LogitN}\\big(\n \\boldsymbol{\\mu} = (${mu[0].toFixed(2)},${mu[1].toFixed(2)}),\n \\boldsymbol{\\sigma} = (${sigma1.toFixed(2)},${sigma2.toFixed(2)}),\n \\rho = ${rho.toFixed(2)}\\big)`}
\n ${logitnormalplt}
\n ${logitNormalColorScale}\n
\n
\n \n","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":716,"value":"

","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":968,"value":"
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n\n---\n#### 🙈 Appendix 🙈\n\nHere is just some data, imports and function defintions that you probably don't want to read.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":972,"value":"#### Plots","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":930,"value":"dirichletplt = ternaryDensity(densityDirichlet, densityResolution, { size: 400, colorScheme: d3.interpolateGreens})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":932,"value":"dirichletColorScale = legend({color: d3.scaleSequential(d3.extent(densityDirichlet.map(d => d.density)), d3.interpolateGreens), title: \"Density\"})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":491,"value":"logitnormalplt = ternaryDensity(densityLogitNormal, densityResolution, { size: 400 })","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":530,"value":"logitNormalColorScale = legend({color: d3.scaleSequential(d3.extent(densityLogitNormal.map(d => d.density)), d3.interpolateBlues), title: \"Density\"})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":758,"value":"#### Imports","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":749,"value":"math = require(\"mathjs\");","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":896,"value":"digamma = require( 'https://cdn.jsdelivr.net/gh/stdlib-js/math-base-special-digamma@umd/browser.js' )","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":898,"value":"trigamma = require( 'https://cdn.jsdelivr.net/gh/stdlib-js/math-base-special-trigamma@umd/browser.js' )","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1036,"value":"import {roundbox, md_colour} from \"@mattiasvillani/imports\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1038,"value":"roundbox","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":762,"value":"#### Data","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":926,"value":"alpha = [inputs.alpha1, inputs.alpha2, inputs.alpha3]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":764,"value":"mu = {\n if (klcloseparam){\n return [digamma(inputs.alpha1) - digamma(inputs.alpha3), digamma(inputs.alpha2) -digamma(inputs.alpha3)]\n }else{\n return [paramlogitn.mu1, paramlogitn.mu2]\n }\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":767,"value":"sigma1 = {\n if (klcloseparam){\n return Math.sqrt(trigamma(inputs.alpha1) + trigamma(inputs.alpha3))\n }else{\n return paramlogitn.sigma1\n }\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":769,"value":"sigma2 = {\n if (klcloseparam){\n return Math.sqrt(trigamma(inputs.alpha2) + trigamma(inputs.alpha3))\n }else{\n return paramlogitn.sigma2\n }\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":771,"value":"rho = {\n if (klcloseparam){\n return trigamma(inputs.alpha3)/(sigma1*sigma2)\n }else{\n return paramlogitn.rho\n }\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1110,"value":"Sigma = {\n if (klcloseparam){\n return [[sigma1**2,rho*sigma1*sigma2],[rho*sigma1*sigma2,sigma2**2]]\n }else{\n return [[paramlogitn.sigma1**2, paramlogitn.rho*paramlogitn.sigma1*paramlogitn.sigma2], \n [paramlogitn.rho*paramlogitn.sigma1*paramlogitn.sigma2, paramlogitn.sigma2**2]]\n }\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":760,"value":"#### Functions","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":747,"value":"function multiLogitNormal_pdf(x, mu, Sigma){\n const p = math.matrix(mu).size()[0]; // Dimension of x is K, dim of mu and Sigma is p = K-1\n \n const logDet = -(p/2)*math.log(2*math.pi) -0.5*math.log(math.det(Sigma));\n const logJacobian = - d3.sum(x.map(d => Math.log(d)))\n //let xtilde = math.multiply( (1/x[p]), x.slice(0, p - 1)) // This is p = K-1 dimensional\n //xtilde = xtilde.map(d => Math.log(d))\n const xtilde = x.slice(0, p).map(d => Math.log(d/x[p]))\n \n const quad = \n -0.5*math.multiply( math.multiply( math.transpose(math.subtract(xtilde,mu)), math.inv(Sigma)), \n math.subtract(xtilde,mu))\n return Math.exp(logDet + logJacobian + quad)\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":919,"value":"import {dirichlet,beta} from \"@mattiasvillani/dirichlet-distribution\"","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":828,"value":"disttext = md`
Ternary density plot from **Herb Susmann** at https://observablehq.com/@herbps10`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":485,"value":"import {ternaryDensity} from \"@mattiasvillani/imports\"","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":914,"value":"densityDirichlet = {\n const density = [];\n for(let i = 1; i < densityResolution; i++) {\n for(let j = 1; j < densityResolution - i; j++) {\n const x1 = i / densityResolution ;\n const x2 = j / densityResolution;\n const x3 = Math.max(0, 1 - x1 - x2);\n const x = [ x1, x2, x3 ];\n\n density.push({ x, density: dirichlet(x, alpha) });\n }\n }\n return density;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":533,"value":"densityLogitNormal = {\n const density = [];\n for(let i = 1; i < densityResolution; i++) {\n for(let j = 1; j < densityResolution - i; j++) {\n const x1 = i / densityResolution ;\n const x2 = j / densityResolution;\n const x3 = Math.max(0, 1 - x1 - x2);\n const x = [ x1, x2, x3 ];\n\n density.push({ x, density: multiLogitNormal_pdf(x, mu, Sigma)});\n }\n }\n return density;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":535,"value":"densityResolution = 40","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":51,"value":"margin = ({ left: 30, top: 30, right: 30, bottom: 30 })","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":8,"value":"height = 400","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":4,"value":"d3 = require(\"d3@v6\", \"d3-hexbin@0.2\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1034,"value":"import {legend} from \"@d3/color-legend\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W33","stdlib":"1","titling":"auto"}} 30:"$2f:body"