2b:null 2c:["$","title",null,{"children":"Bayesian credible intervals | Mattias Villani | Observable"}] 31:T64b,A 95% Bayesian posterior **credible interval** ${tex`[l,u]`} for a parameter ${tex`\theta`} is an interval that contains 95% of the posterior probability mass *conditional on the observed data* ${tex`\boldsymbol{x}`}: ${tex.block`\mathrm{Pr}(\theta \in [l,u] \mid \boldsymbol{x}) = 0.95`} Note that the interpretation of a credible interval is different from a frequentist confidence interval: - a credible interval is really the (posterior) probability that the unknown parameter ${tex`\theta`} is in the interval ${tex`[l,u]`} 😍. This probability is *conditional on the single data set ${tex`\boldsymbol{x}`} that we have actually observed*. - a frequentist [confidence interval](confidence-interval-for-a-proportion) is a random interval that covers the population parameter ${tex`\theta`} *in 95% of all datasets* drawn from the population 😵‍💫. There are many ways to construct a credibility interval with a certain coverage probability. This widget compares two types of intervals: - An **equal tail** interval cuts off equal probability in the left and right tail; for example, a 95% interval sets *l* and *u* to the 2.5% and 97.5% posterior quantiles, respectively. - A 95% Highest Posterior Density (HPD) interval is made up of the 95% most probable parameter values in the posterior. HPD intervals can be disconnected, so **HPD regions** is a more appropriate name. There is of course nothing special about 95% coverage, the coverage of the interval can be anything between 0% and 100%. The widget below allows you to experiment with different distributions and coverage probabilities. 32:Ta97,viewof params = { var form; if (dist_type == 'beta'){ form = Inputs.form([ Inputs.range([0.1, 20], {value: 3, step: 0.1, label: tex`\text{shape }\alpha`}), Inputs.range([0.1, 20], {value: 5, step: 0.1, label: tex`\text{shape }\beta`}) ]) } else if (dist_type == 'exponential'){ form = Inputs.form([ Inputs.range([0.01, 10], {value: 1, step: 0.01, label: tex`\text{rate }\lambda`}) ]) } else if (dist_type == 'studentt'){ form = Inputs.form([ Inputs.range([-5, 5], {value: 0, step: 0.1, label: tex`\text{location }\mu`}), Inputs.range([0.1, 5], {value: 1, step: 0.1, label: tex`\text{scale }\sigma`}), Inputs.range([1, 50], {value: 4, step: 1, label: tex`\text{degrees of freedom }\nu`}) ]) } else if (dist_type == 'chi2'){ form = Inputs.form([ Inputs.range([1, 50], {value: 3, step: 1, label: tex`\text{degrees of freedom }\nu`}) ]) } else if (dist_type == 'lognormal'){ form = Inputs.form([ Inputs.range([-5,5], {value: 0, step: 0.1, label: tex`\text{location }\mu`}), Inputs.range([0.1,5], {value: 1, step: 0.1, label: tex`\text{scale }\sigma`}) ]) } else if (dist_type == 'normal'){ form = Inputs.form([ Inputs.range([-5,5], {value: 0, step: 0.1, label: tex`\text{mean }\mu`}), Inputs.range([0.1,5], {value: 1, step: 0.1, label: tex`\text{stdev }\sigma`}) ]) } else if (dist_type == 'binormal mixture'){ form = Inputs.form([ Inputs.range([-5,5], {value: -2, step: 0.1, label: tex`\text{location }\mu_1`}), Inputs.range([-5,5], {value: 2, step: 0.1, label: tex`\text{location }\mu_2`}), Inputs.range([0.1,5], {value: 1, step: 0.1, label: tex`\text{scale }\sigma_1`}), Inputs.range([0.1,5], {value: 1, step: 0.1, label: tex`\text{scale }\sigma_2`}), Inputs.range([0,1], {value: 0.3, step: 0.01, label: tex`\text{weight }w_1`}) ]) } else if (dist_type == 'trinormal mixture'){ form = Inputs.form([ Inputs.range([-5,5], {value: -2, step: 0.1, label: tex`\text{location }\mu_1`}), Inputs.range([-5,5], {value: 0, step: 0.1, label: tex`\text{location }\mu_2`}), Inputs.range([-5,5], {value: 2, step: 0.1, label: tex`\text{location }\mu_3`}), Inputs.range([0.1,5], {value: 0.5, step: 0.1, label: tex`\text{scale }\sigma_1`}), Inputs.range([0.1,5], {value: 0.5, step: 0.1, label: tex`\text{scale }\sigma_2`}), Inputs.range([0.1,5], {value: 0.5, step: 0.1, label: tex`\text{scale }\sigma_3`}), Inputs.range([0,0.49], {value: 0.3, step: 0.01, label: tex`\text{weight }w_1`}), Inputs.range([0,0.5], {value: 0.3, step: 0.01, label: tex`\text{weight }w_2`}) ]) } return form; }33:T805,distribution = { let d; if (dist_type == "normal"){ d = ({ "pdf" : x => jStat.normal.pdf(x,params[0],params[1]), "quantile" : x => jStat.normal.inv(x,params[0],params[1]), "support" : [-Infinity,Infinity] }) } else if (dist_type == "binormal mixture"){ d = ({ "pdf" : x => params[4]*jStat.normal.pdf(x,params[0],params[2]) + (1-params[4])*jStat.normal.pdf(x,params[1],params[3]), "quantile" : x => mon2inv(x), "support" : [-Infinity,Infinity] }) } else if (dist_type == "trinormal mixture"){ d = ({ "pdf" : x => params[6]*jStat.normal.pdf(x,params[0],params[3]) + params[7]*jStat.normal.pdf(x,params[1],params[4]) + (1 -params[6] - params[7])*jStat.normal.pdf(x,params[2],params[5]), "quantile" : x => mon3inv(x), "support" : [-Infinity,Infinity] }) } else if (dist_type == "chi2"){ d = ({ "pdf" : x => jStat.chisquare.pdf(x,params[0]), "quantile" : x => jStat.chisquare.inv(x,params[0]), "support" : [Number.EPSILON,Infinity] }) } else if (dist_type == "exponential"){ d = ({ "pdf" : x => jStat.exponential.pdf(x,params[0]), "quantile" : x => jStat.exponential.inv(x,params[0]), "support" : [Number.EPSILON,Infinity] }) } else if (dist_type == "lognormal"){ d = ({ "pdf" : x => jStat.lognormal.pdf(x,params[0],params[1]), "quantile" : x => jStat.lognormal.inv(x,params[0],params[1]), "support" : [Number.EPSILON,Infinity] }) } else if (dist_type == "beta"){ d = ({ "pdf" : x => jStat.beta.pdf(x,params[0],params[1]), "quantile" : x => jStat.beta.inv(x,params[0],params[1]), "support" : [Number.EPSILON,1-Number.EPSILON] }) } else if (dist_type == "studentt"){ d = ({ "pdf" : x => jStat.studentt.pdf( (x-params[0])/params[1], params[2])/params[1], "quantile" : x => params[0] + params[1]*jStat.studentt.inv(x,params[2]), "support" : [-Infinity,Infinity] }) } return d }34:T728,function HPDregions(distribution, coverage) { const quantiles = [0.001, 0.999]; const minMax = quantiles.map(quantile => distribution.quantile(quantile)); // Check if end of support has higher density const minDensity = distribution.pdf(minMax[0]); const supportMin = distribution.support[0]; if (minDensity < distribution.pdf(supportMin)) { minMax[0] = supportMin; } const maxDensity = distribution.pdf(minMax[1]); const supportMax = distribution.support[1]; if (maxDensity < distribution.pdf(supportMax)) { minMax[1] = supportMax; } const xGrid = d3.range(minMax[0], minMax[1], (minMax[1] - minMax[0])/1000); const binSize = xGrid[1] - xGrid[0]; const dens = xGrid.map(x => distribution.pdf(x)); const sortedIndices = dens.map((_, i) => i).sort((a, b) => dens[b] - dens[a]); const xSort = sortedIndices.map(i => xGrid[i]); const densSort = sortedIndices.map(i => dens[i]); let finalPointIn; let cumulativeSum = 0; for (let i = 0; i < densSort.length; i++) { cumulativeSum += densSort[i] * binSize; if (cumulativeSum >= coverage) { finalPointIn = i; break; } } const actualCoverage = cumulativeSum; const hpdPoints = d3.sort(xSort.slice(0, finalPointIn + 1)); const breakpoints = [0]; for (let i = 1; i < hpdPoints.length; i++) { if (hpdPoints[i] - hpdPoints[i - 1] > 1.9 * binSize) { breakpoints.push(i-1); } } breakpoints.push(hpdPoints.length - 1 ); const nIntervals = breakpoints.length - 1; const hpd = new Array(nIntervals); for (let j = 0; j < nIntervals; j++) { hpd[j] = [hpdPoints[breakpoints[j] + 1], hpdPoints[breakpoints[j + 1]]]; } return [hpd, actualCoverage, hpdPoints] } 35:T49d,hpdtext = { const numprec = 3; switch (hpd[0].length){ case 1: return(md`The ${coverage}% HPD region is (${hpd[0][0][0].toFixed(numprec)}, ${hpd[0][0][1].toFixed(numprec)})
The ${coverage}% equal tail region is (${equaltailinterval[0].toFixed(numprec)}, ${equaltailinterval[1].toFixed(numprec)})`) case 2: return(md`The ${coverage}% HPD region is (${hpd[0][0][0].toFixed(numprec)}, ${hpd[0][0][1].toFixed(numprec)}) ${tex`\cup`} (${hpd[0][1][0].toFixed(numprec)}, ${hpd[0][1][1].toFixed(numprec)})
The ${coverage}% equal tail region is (${equaltailinterval[0].toFixed(numprec)}, ${equaltailinterval[1].toFixed(numprec)})`) case 3: return(md`The ${coverage}% HPD region is (${hpd[0][0][0].toFixed(numprec)}, ${hpd[0][0][1].toFixed(numprec)}) ${tex`\cup`} (${hpd[0][1][0].toFixed(numprec)}, ${hpd[0][1][1].toFixed(numprec)}) ${tex`\cup`} (${hpd[0][2][0].toFixed(numprec)}, ${hpd[0][2][1].toFixed(numprec)})
The ${coverage}% equal tail region is (${equaltailinterval[0].toFixed(numprec)}, ${equaltailinterval[1].toFixed(numprec)})`) default: return("hee") } }36:[] 2f:{"id":"6d25175cf8e61aab","slug":"bayesian-credible-intervals","title":"Bayesian credible intervals","thumbnail":"34b7cc3e8976fab37cb9ab9c1b4ae7356abe7c2bf75c7a3e64fd038592a55d56","access_level":"public","create_time":"$D2024-02-03T10:25:15.615Z","update_time":"$D2025-09-11T12:43:24.046Z","view_time":null,"first_public_time":"$D2024-02-03T16:25:35.832Z","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":596,"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":1,"collection":{"id":"9fdb7c9ba7be6ba8","slug":"bayesian-learning","title":"Bayesian Learning","owner":{"id":"7592a7ab121f74a0","login":"mattiasvillani","name":"Mattias Villani"}},"type":"notebook","body":{"title":"Bayesian credible intervals","cells":[{"id":0,"value":"# Bayesian credible intervals","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":514,"value":"$31","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":582,"value":"anchor = md`### Interactive exploration of credible intervals`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":23,"value":"viewof dist_type = Inputs.select(['normal', 'binormal mixture', 'trinormal mixture', 'beta', 'chi2', 'exponential', 'lognormal','studentt'], {value: \"binormal\", label: \"distribution:\"})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":26,"value":"$32","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":163,"value":"viewof coverage = Inputs.range([1, 99], {value: 90, label: \"coverage %\", step: 1})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":378,"value":"$${hpdtext}","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":124,"value":"Plot.plot({\ncolor: {\n legend: true,\n domain: [\"posterior density\",\"HPD region\", \"equal tail interval\"],\n range: [colors[2], d3.rgb(192, 163, 77, 0.3), colors[0]]\n },\n y: {\n label: \"density\"\n },\n x: {\n label: \"x\"\n },\n marks: [\n //Plot.ruleX([d3.min(pdfdata.map(d => d.x))]),\n Plot.areaY(pdfdata, {filter: d => HPDregion_number(d.x) == 1, x: \"x\", y: \"pdf\", fill: colors[1], opacity: 0.3}),\n Plot.areaY(pdfdata, {filter: d => HPDregion_number(d.x) == 2, x: \"x\", y: \"pdf\", fill: colors[1], opacity: 0.3}),\n Plot.areaY(pdfdata, {filter: d => HPDregion_number(d.x) == 3, x: \"x\", y: \"pdf\", fill: colors[1], opacity: 0.3}),\n Plot.areaY(pdfdata, {filter: d => HPDregion_number(d.x) == 4, x: \"x\", y: \"pdf\", fill: colors[1], opacity: 0.3}),\n Plot.lineY(pdfdata, {x: \"x\", y: \"pdf\", stroke: colors[2]}),\n Plot.line([{x: equaltailinterval[0], y:0}, {x: equaltailinterval[1], y: 0}],{x: \"x\", y: \"y\", stroke: colors[0], strokeWidth: 4})\n ]\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":5,"value":"
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n\n---\n## 🙈 Appendix 🙈\n\nJust some code and implementation details that you probably don't want to read.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":18,"value":"$33","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":11,"value":"#### Imports","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":15,"value":"jStat = require('jstat')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":465,"value":"import {NewtonsMethod, secant} from \"@mattiasvillani/statistics-tools\"","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":7,"value":"#### Functions","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":9,"value":"$34","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":139,"value":"function HPDregion_number(x){\n let region_number = 0;\n let count = 0;\n for (const hpdregion of hpd[0]) {\n count = count + 1;\n if (x >= hpdregion[0] && x <= hpdregion[1]){\n region_number = count;\n break;\n }\n }\n return region_number\n}\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":293,"value":"function approxquant2(p){\n if (p < 0.5){\n return d3.min([jStat.normal.inv(p, params[0],params[2]),jStat.normal.inv(p, params[1],params[3])])\n }\n else{\n return d3.max([jStat.normal.inv(p, params[0],params[2]),jStat.normal.inv(p, params[1],params[3])])\n }\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":320,"value":"function approxquant3(p){\n if (p < 0.5){\n return d3.min([jStat.normal.inv(p, params[0],params[3]),\n jStat.normal.inv(p, params[1],params[4]),\n jStat.normal.inv(p, params[2],params[5])])\n }\n else{\n return d3.max([jStat.normal.inv(p, params[0],params[3]),\n jStat.normal.inv(p, params[1],params[4]),\n jStat.normal.inv(p, params[2],params[5])])\n }\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":119,"value":"#### Data","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":326,"value":"colors = [\"#6C8EBF\", \"#c0a34d\", \"#780000\", \"#007878\", \"#bf9d6c\", \"#3A6B35\", \"#b5c6df\",\"#eadaaa\", \"#bb989a\", \"#98bbb9\", \"#bf8d6c\", \"#CBD18F\"]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":185,"value":"hpd = HPDregions(distribution, coverage/100)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":351,"value":"$35","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":121,"value":"pdfdata = {\n const lower = isFinite(distribution.support[0])? distribution.support[0] : distribution.quantile(0.001)\n const upper = isFinite(distribution.support[1])? distribution.support[1] : distribution.quantile(0.999)\n const x = d3.range(lower, upper, (upper-lower)/10000)\n\n const pdf = x.map(x => ({x: x, pdf: distribution.pdf(x), hpdregion: HPDregion_number(x)}));\n return pdf\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":389,"value":"equaltailinterval = [distribution.quantile((1-coverage/100)/2), distribution.quantile((1+coverage/100)/2)]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":458,"value":"function mon_cdf2(x, p){\n return params[4]*jStat.normal.cdf(x,params[0],params[2]) + (1-params[4])*jStat.normal.cdf(x,params[1],params[3]) - p\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":468,"value":"function mon_pdf2(x){\n return params[4]*jStat.normal.pdf(x,params[0],params[2]) + (1-params[4])*jStat.normal.pdf(x,params[1],params[3])\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":461,"value":"function mon_cdf3(x,p){\n return params[6]*jStat.normal.cdf(x,params[0],params[3]) + \n params[7]*jStat.normal.cdf(x,params[1],params[4]) + \n (1 -params[6] - params[7])*jStat.normal.cdf(x,params[2],params[5]) - p\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":470,"value":"function mon_pdf3(x){\n return params[6]*jStat.normal.pdf(x,params[0],params[3]) + \n params[7]*jStat.normal.pdf(x,params[1],params[4]) + \n (1 -params[6] - params[7])*jStat.normal.pdf(x,params[2],params[5])\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":473,"value":"function mon2inv(p){\n return NewtonsMethod(x => mon_cdf2(x, p), mon_pdf2, approxquant2(p))\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":479,"value":"function mon3inv(p){\n return NewtonsMethod(x => mon_cdf3(x, p), mon_pdf3, approxquant3(p))\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W36","stdlib":"1","titling":"auto"}} 30:"$2f:body"