2b:null 2c:["$","title",null,{"children":"Posterior approximation - Beta model for proportions | Mattias Villani | Observable"}] 31:T439,Consider data in the form of proportions following an iid Beta distribution ${tex`X_1,\ldots,X_n \overset{\mathrm{iid}}{\sim} \mathrm{Beta}(\alpha,\beta)`} The likelihood function is ${tex` p(\boldsymbol{x} \vert \alpha,\beta) = \prod_{i=1}^n \frac{1}{\mathrm{B}(\alpha,\beta)} x_i^{\alpha -1}(1-x_i)^{\beta -1} = \frac{1}{\mathrm{B(\alpha,\beta)}^n} (\prod_{i=1}^n x_i)^{\alpha -1}(\prod_{i=1}^n (1-x_i))^{\beta -1} `} and the log-likelihood function is ${tex` \ell(\alpha,\beta)=\log p(\boldsymbol{x} \vert \alpha,\beta) = -n \log B(\alpha,\beta) + (\alpha - 1)(\sum_{i=1}^n \log x_i) + (\beta -1)(\sum_{i=1}^n \log(1-x_i)) `} We will here do a normal approximation of the posterior distribution ${tex`\alpha,\beta \vert x_1,\ldots,x_n \sim N(\boldsymbol{\mu}, \boldsymbol{\Sigma})`}, where ${tex`\boldsymbol{\mu} = (\tilde\alpha,\tilde\beta)^\top`} is the posterior mode and ${tex`\boldsymbol{\Sigma} =-\boldsymbol{H}^{-1}`}, where ${tex`\boldsymbol{H}`} is the Hessian matrix evaluated at the mode. We assume a flat prior ${tex`p(\alpha,\beta)\propto \text{constant}`}.32:[] 2f:{"id":"aa335de8b718294a","slug":"posterior-approximation-beta-model","title":"Posterior approximation - Beta model for proportions","thumbnail":"05f33174e952c07f598503177b606efd35cd09f5f302908ab31d967a2c13ca57","access_level":"public","create_time":"$D2023-11-18T22:48:13.965Z","update_time":"$D2025-02-11T17:16:05.644Z","view_time":null,"first_public_time":"$D2023-11-19T22:21:01.816Z","fork":false,"stars":0,"starred":false,"owner":{"type":"team","id":"7592a7ab121f74a0","login":"mattiasvillani","name":"Mattias Villani","avatar_url":"https://avatars-next.observableusercontent.com/a440aaf39a3abefb86cbcc758fd7a5ab42555ae68d6825b52ec006f385a94d1d","tier":"pro_2024"},"version":563,"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":true,"trashed":false,"subscribed":null,"forks":1,"collection":{"id":"9fdb7c9ba7be6ba8","slug":"bayesian-learning","title":"Bayesian Learning","owner":{"id":"7592a7ab121f74a0","login":"mattiasvillani","name":"Mattias Villani"}},"type":"notebook","body":{"title":"Posterior approximation - Beta model for proportions","cells":[{"id":0,"value":"# Posterior approximation - Beta model for proportions","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":3,"value":"$31","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":399,"value":"\n**Gradient** of log-likelihood\n\n${tex`\n \\begin{pmatrix}\n \\frac{\\partial \\log p(\\boldsymbol{x} \\vert \\alpha,\\beta)}{\\partial \\alpha} \\\\\n \\frac{\\partial \\log p(\\boldsymbol{x} \\vert \\alpha,\\beta)}{\\partial \\beta}\n \\end{pmatrix}\n = \n \\begin{pmatrix}\n -n\\big(\\psi^{(0)}(\\alpha) - \\psi^{(0)}(\\alpha + \\beta)\\big) + \\sum_{i=1}^n \\log x_i \\\\\n -n\\big(\\psi^{(0)}(\\beta) - \\psi^{(0)}(\\alpha + \\beta)\\big) + \\sum_{i=1}^n \\log (1-x_i)\n \\end{pmatrix}\n`}\n\nwhere ${tex`\\psi^{(k)}(z):=\\frac{\\partial^{k+1}}{\\partial z^{k+1}}\\log\\Gamma(z)`} is the [polygamma](https://en.wikipedia.org/wiki/Polygamma_function) function of order *k*.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":185,"value":"**Hessian** matrix\n\n${tex`\n\\boldsymbol{H} = \n\\begin{pmatrix} \n \\frac{\\partial^2}{\\partial\\alpha^2}\\ell(\\alpha,\\beta) & \\frac{\\partial}{\\partial\\alpha\\partial\\beta}\\ell(\\alpha,\\beta)\\\\ \n \\frac{\\partial}{\\partial\\alpha\\partial\\beta}\\ell(\\alpha,\\beta) & \\frac{\\partial^2}{\\partial\\beta^2}\\ell(\\alpha,\\beta)\n\\end{pmatrix}\n=\n-n\\begin{pmatrix} \n \\psi^{(1)}(\\alpha) - \\psi^{(1)}(\\alpha + \\beta) & -\\psi^{(1)}(\\alpha + \\beta)\\\\ \n -\\psi^{(1)}(\\alpha + \\beta) & \\psi^{(1)}(\\beta) - \\psi^{(1)}(\\alpha + \\beta) \n\\end{pmatrix}\n\n`}","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":381,"value":"anchor = md`The interactive graph below gives the contours of the multivariate normal approximation for data simulated from the model.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":512,"value":"viewof n = Inputs.range([1, max_sample_size], {value: 50, label: \"Sample size\", step: 1})\n ","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":30,"value":"viewof dgp_settings = Inputs.form([\n Inputs.range([0.1, 30], {value: 10, label: tex`\\alpha`, step: 0.1}),\n Inputs.range([0.1, 30], {value: 10, label: tex`\\beta`, step: 0.1})\n ])","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":488,"value":"\n \n \n \n

\n ${plttrue}\n

\n ${pltapprox}\n

\n \n","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":325,"value":"
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n\n---\n## 🙈 Appendix 🙈\nHere follow code and details that you probably do not want to read.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":423,"value":"plttrue = Plot.plot({\n x: {label: \"α\", domain: [0,mle[0] + 5*stdev_alpha]},\n y: {label: \"β\", domain: [0,mle[1] + 5*stdev_beta]},\n color: {\n legend: true,\n type: \"diverging\",\n label: \"true posterior density\",\n //domain: [0,1]\n },\n marks: [\n Plot.contour(loglik_on_grid, {filter: d => d.lik > Number.EPSILON, x: \"alpha\", y: \"beta\", fill: \"lik\", stroke: \"black\", blur: 2})\n ]\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":308,"value":"pltapprox = Plot.plot({\n x: {label: \"α\", domain: [0,mle[0] + 5*stdev_alpha]},\n y: {label: \"β\", domain: [0,mle[1] + 5*stdev_beta]},\n color: {\n //domain: [0, 0.01],\n legend: true,\n type: \"diverging\",\n label: \"density of normal posterior approximation\"\n },\n marks: [\n Plot.contour(mvnormal_on_grid, {filter: d => d.postpdfapprox > Number.EPSILON, x: \"alpha\", y: \"beta\", fill: \"postpdfapprox\", stroke: \"black\", blur: 2})\n ]\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":326,"value":"#### Data","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":334,"value":"max_sample_size = 100","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":49,"value":"alpha = dgp_settings[0]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":51,"value":"beta = dgp_settings[1]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":42,"value":"data = d3.range(max_sample_size).map(d => jstat.beta.sample(alpha, beta))","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":329,"value":"x = data.slice(0,n)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":36,"value":"sumlogx = d3.sum(x.map(d => Math.log(d)))","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":40,"value":"sumlog1minusx = d3.sum(x.map(d => Math.log(1-d)))","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":16,"value":"function loglikfunc(param){\n let alpha = param[0]\n let beta = param[1]\n return -n*jstat.betaln( alpha, beta ) + (alpha-1)*sumlogx + (beta-1)*sumlog1minusx\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":431,"value":"function negloglikfunc(param){\n let alpha = param[0]\n let beta = param[1]\n return -(-n*jstat.betaln( alpha, beta ) + (alpha-1)*sumlogx + (beta-1)*sumlog1minusx)\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":65,"value":"mle_opt = optimjs.minimize_Powell(negloglikfunc, [16,12])","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":105,"value":"mle = mle_opt.argument","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":109,"value":"mle_max = mle_opt.fncvalue","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":133,"value":"function neghessian(param){\n const alpha = param[0]\n const beta = param[1]\n const h11 = n*(trigamma(alpha) - trigamma(alpha + beta))\n const h22 = n*(trigamma(beta) - trigamma(alpha + beta))\n const h12 = -n*trigamma(alpha+beta)\n return [[h11,h12],[h12,h22]]\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":138,"value":"Sigma = math.inv(neghessian([alpha,beta]))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":279,"value":"stdev_alpha = Math.sqrt(Sigma[0][0])","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":283,"value":"stdev_beta = Math.sqrt(Sigma[1][1])","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":163,"value":"maxpdf = mvnpdf([mle[0],mle[1]], [mle[0],mle[1]], Sigma)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":216,"value":"loglik_on_grid = {\n const ngridalpha = 30;\n const ngridbeta = 30;\n var griddata = [];\n for (let alpha of d3.range(d3.max([0, mle[0] - 3*stdev_alpha]), mle[0] + 3*stdev_alpha, 6*stdev_alpha/ngridalpha)){\n for (let beta of d3.range(d3.max([0, mle[1] - 3*stdev_beta]), mle[1] + 3*stdev_beta, 6*stdev_beta/ngridbeta)){\n griddata.push({alpha: alpha, beta: beta, loglik: loglikfunc([alpha,beta])})\n //griddata.push({alpha: alpha, beta: beta, loglik: alpha+beta})\n }\n }\n const binarea = (6*stdev_alpha/ngridalpha)*(6*stdev_beta/ngridbeta);\n const max_on_grid = d3.max(griddata.map(d => d.loglik));\n const normconst = binarea*d3.sum(griddata.map(d => Math.exp(d.loglik - max_on_grid)))\n return griddata.map(d => ({alpha: d.alpha, beta: d.beta, lik: Math.exp(d.loglik - max_on_grid)/normconst}))\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":265,"value":"mvnormal_on_grid = {\n const ngridalpha = 30;\n const ngridbeta = 30;\n var griddata = [];\n for (let alpha of d3.range(d3.max([0, mle[0] - 3*stdev_alpha]), mle[0] + 3*stdev_alpha, 6*stdev_alpha/ngridalpha)){\n for (let beta of d3.range(d3.max([0, mle[1] - 3*stdev_beta]), mle[1] + 3*stdev_beta, 6*stdev_beta/ngridbeta)){\n griddata.push({alpha: alpha, beta: beta, postpdfapprox: mvnpdf([alpha, beta], [mle[0],mle[1]], Sigma) })\n //griddata.push({alpha: alpha, beta: beta, loglik: alpha+beta})\n }\n }\n return griddata\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":27,"value":"##### Imports","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":25,"value":"jstat = require('jstat')","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":135,"value":"math = require(\"mathjs\");","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":69,"value":"optimjs = require('optimization-js@1.5.0/dist/optimization.js')","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":21,"value":"digamma = require( 'https://cdn.jsdelivr.net/gh/stdlib-js/math-base-special-digamma@umd/browser.js' )","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":23,"value":"trigamma = require( 'https://cdn.jsdelivr.net/gh/stdlib-js/math-base-special-trigamma@umd/browser.js' )","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":155,"value":"import {mvnpdf} from \"@mattiasvillani/multivariate-normal-distribution\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W32","stdlib":"1","titling":"auto"}} 30:"$2f:body"