29:null 2a:["$","title",null,{"children":"Pólya-Gamma distribution | Mattias Villani | Observable"}] 2f:[] 2d:{"id":"de60acdc96880667","slug":"polya-gamma-distribution","title":"Pólya-Gamma distribution","thumbnail":"374dbdb719017ea61ffbf68793ab2330c82b54885124e18caf8c66719f916ecf","access_level":"public","create_time":"$D2023-11-07T20:37:37.524Z","update_time":"$D2023-11-08T09:05:37.204Z","view_time":null,"first_public_time":"$D2023-11-07T21:44:11.933Z","has_unpublished_edits":null,"fork":{"id":"811f9749bca2faf5","slug":"skew-normal-distribution","title":"Skew-normal distribution","owner":{"id":"7592a7ab121f74a0","login":"mattiasvillani","name":"Mattias Villani"}},"stars":0,"starred":false,"owner":{"type":"team","id":"7592a7ab121f74a0","login":"mattiasvillani","name":"Mattias Villani","avatar_url":"https://avatars-next.observableusercontent.com/a440aaf39a3abefb86cbcc758fd7a5ab42555ae68d6825b52ec006f385a94d1d","tier":"education_2026"},"version":1107,"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,"lapsed":null,"subscribed":null,"forks":0,"collection":{"id":"40c8ef5af1a51aff","slug":"distributions","title":"Statistical distributions","owner":{"id":"7592a7ab121f74a0","login":"mattiasvillani","name":"Mattias Villani"}},"type":"notebook","body":{"title":"Pólya-Gamma distribution","cells":[{"id":38,"value":"# Pólya-Gamma distribution","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":629,"value":"

${tex`X\\sim \\mathrm{PG}(b,c)`}, for ${tex`X>0`}","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":496,"value":"A ${tex`X\\sim \\mathrm{PG}(b,c)`} random variable is defined as a convolution of Gamma variables\n\n${tex`\\qquad X \\overset{d}{=} \\sum_{k=1}^\\infty v_k Y_k`}\n\nwhere ${tex`\\overset{d}{=}`} mean equality in distribution, ${tex`Y_k \\overset{iid}{\\sim} \\mathrm{Gamma}(b,1)`} and the weights are \n\n${tex`\\qquad v_k = \\frac{1}{2(k-1/2)^2\\pi^2 +c^2/2}`}.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":118,"value":"viewof inputs = Inputs.form([\n Inputs.range([1, 10], {value: 1, step: 1, label: tex`b`}),\n Inputs.range([0, 10], {value: 1, step: 0.1, label: tex`c`}),\n Inputs.range([0, 5], {value: 0.2, step: 0.01, label: \"quantile:\"})\n ])","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":461,"value":"dist_quantile = tex`\\text{If } X \\sim \\mathrm{PG}(${b},${c.toPrecision(2)}) \\text{ then } `","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":555,"value":"moments = {\n return tex`\n \\begin{aligned}\n &E( X) = \\frac{b}{2c}\\Big( \\frac{e^c -1}{1 + e^c} \\Big) = ${mean.toPrecision(3)} \\\\[1em] \n &\\mathrm{Stdev}( X) = \\sqrt{ (b / 4c^3)(\\sinh(c) - c)(\\mathrm{sech}(c/2)^2) }= \n ${(Math.sqrt(variance)).toPrecision(3)}\\\\[1em]\n &P( X\\leq ${quantile})=${cdf.toPrecision(4)}\n \\end{aligned}\n `\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":305,"value":"plt_pdf = Plot.plot({\n color: {\n legend: true\n },\n x: {\n label: \"x\",\n axis: true,\n },\n y: {\n label: \"f(x)\",\n //axis: false,\n //domain: [0,jstat.chisquare.pdf(mode, nu1, nu2)]\n },\n marks: [\n Plot.ruleY([0]),\n Plot.ruleX([0]),\n Plot.line(pdfdata, {x: \"x\", y: \"pdf\", stroke : \"purple\", strokeWidth: 2}),\n Plot.areaY(pdfdata, {filter: d => d.x <= quantile, x: \"x\", y: \"pdf\", fill: \"purple\", opacity: 0.2})\n ]\n })","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":109,"value":"
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n\n---\n## 🙈 Appendix 🙈\n\nHere is some code and data details that you probably do not want to read.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":221,"value":"#### Imports","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":26,"value":"jstat = require('jstat')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":750,"value":"mathjs = require('mathjs')","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":223,"value":"#### Data","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":965,"value":"b = inputs[0]","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":966,"value":"c = inputs[1]","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":662,"value":"quantile = inputs[2]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":706,"value":"mean = (b/(2*c))*((Math.exp(c) - 1 )/(1 + Math.exp(c)))","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":766,"value":"variance = (b / (4 * c**3)) * (Math.sinh(c) - c) * (mathjs.sech(c/2)**2)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":897,"value":"gridsize = 0.001","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1019,"value":"TRUNC = 0.64","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1022,"value":"cutoff = 1 / TRUNC","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1024,"value":"ntrunc = 20","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1052,"value":"pgpdf(0.4, 1, 1)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":830,"value":"function pgpdf(x, b, c){\n if (b == 1){\n return Math.exp(jacobi_logpdf(c/2, 4*x, ntrunc) + Math.log(4))\n }\n else{\n return Math.exp(pg_logpdf(b, c, x, ntrunc))\n }\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":974,"value":"function jacobi_logpdf(z, x, ntrunc){\n var v = 0;\n for (let n = 0; n <= ntrunc ; n++){\n v += (n % 2 === 0 ? 1 : -1) * acoef(n, x)\n }\n return logcosh(z) -x*(z**2)/2 + Math.log(v)\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":987,"value":"function acoef(n, x, r = TRUNC){\n if (x>TRUNC){\n return Math.PI * (n + 0.5) * Math.exp( -((n+0.5)**2) * Math.PI**2 * x/2)\n }else{\n return (2/Math.PI/x)**1.5 * Math.PI * (n + 0.5) * Math.exp( -2*(n+0.5)**2/x )\n }\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":991,"value":"function logcosh(x){\n return x + Math.log1p(Math.exp(-2*x)) - Math.log(2) \n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1030,"value":"function pg_logpdf(b, c, x, ntrunc){\n return b*logcosh(c/2) - x*c**2/2 + pg0_logpdf(x, b, ntrunc)\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1032,"value":"function pg0_logpdf(x, b, ntrunc){\n var v = 0;\n for (let n = 0; n <= ntrunc ; n++){\n v += (n % 2 === 0 ? 1 : -1) * Math.exp(pg_logcoef(x, b, n))\n }\n return (b-1)*Math.log(2) - mathjs.lgamma(b) + Math.log(v)\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1036,"value":"function pg_logcoef(x, b, n){\n return mathjs.lgamma(n+b)-mathjs.lgamma(n+1) + Math.log(2*n+b) - Math.log(2*Math.PI * x**3)/2 - ((2*n+b)**2)/(8*x)\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":17,"value":"pdfdata = {\n const x = d3.range(Number.EPSILON, 3, gridsize); \n const pdf = x.map(x => ({x: x, pdf: pgpdf(x, b, c)}));\n return pdf\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":19,"value":"cdf = d3.sum(pdfdata.filter(d => d.x <=quantile).map(d => d.pdf))*gridsize","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":229,"value":"#### Style","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":202,"value":"","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":58,"value":"viewof md_colour = Inputs.color({label: \"Highlight color cells\", value: \"#F5F5DC\"})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W2f","stdlib":"1","titling":"auto"}} 2e:"$2d:body"