2b:null 2c:["$","title",null,{"children":"Coverage of interval estimates for a proportion | Mattias Villani | Observable"}] 32:T572,Plot.plot({ caption: plotValidProb ? html` The light gray area marks out the p where np<10 or n(1-p)<10
The light purple line plots the probability that np̂≥10 and n(1-p̂)≥10 for the estimator p̂=s/n. ` : html``, color: { legend: true, domain: ["wald", "wilson", "bayes"], range: [colors[0], colors[1], colors[2]] }, y: { label: "coverage", domain: [lower_y,1] }, x: { label: "p", domain: [0,1] }, marks: [ // Plot.ruleX(coveragedata, {x: "p", y: "conflevel", stroke: "gray", strokeWidth: 0.5}), Plot.areaY(coveragedata, {filter: plotValidProb, x: "p", y: "invalidp", fill: "lightgray", opacity: 0.3}), Plot.ruleY([confidence/100], {stroke: "gray", strokeWidth: 1}), Plot.ruleY([0], {stroke: "black", strokeWidth: 1}), Plot.ruleX([0], {stroke: "black", strokeWidth: 1}), Plot.lineY(coveragedata, {filter: plotValidProb, x: "p", y: "prValid", stroke: "#800080", strokeWidth: 2}), Plot.lineY(coveragedata, {filter: plotselector.includes("wald"), x: "p", y: "wald", stroke: colors[0]}), Plot.lineY(coveragedata, {filter: plotselector.includes("wilson"), x: "p", y: "wilson", stroke: colors[1]}), Plot.lineY(coveragedata, {filter: plotselector.includes("bayes"), x: "p", y: "bayes", stroke: colors[2]}), ] })33:T464,coveragedata = { let coverages = [] const sgrid = d3.range(0, n+1) for (let i = 0; i < pgrid.length; ++i) { let p = pgrid[i] let invalidp = (n*p < 10) || (n*(1-p) < 10) ? 1 : 0; let prValid = (n >= 20) ? jstat.binomial.cdf( n-10, n, p ) - jstat.binomial.cdf( 9, n, p ) : 0 let binomprob = sgrid.map(s => jstat.binomial.pdf(s, n, p)) let coverageWald = 0; let coverageWilson = 0; let coverageBayes = 0; for (let s = 0; s <= n; ++s) { let waldinterval = wald(s, n, zcrit) if (waldinterval[0] <= p && waldinterval[1] >= p){coverageWald += binomprob[s]} let wilsoninterval = wilson(s, n, zcrit) if (wilsoninterval[0] <= p && wilsoninterval[1] >= p){coverageWilson += binomprob[s]} let bayesinterval = bayesian(s, n, confidence, alpha, beta) if (bayesinterval[0] <= p && bayesinterval[1] >= p){coverageBayes += binomprob[s]} } coverages.push(({p: p, invalidp: invalidp, prValid: prValid, wald: coverageWald, wilson: coverageWilson, bayes: coverageBayes, conflevel: confidence/100, priorpdf: jstat.beta.pdf(p,alpha,beta)})) } return coverages }34:[] 2f:{"id":"a7b18022453fe4ef","slug":"coverage-of-confidence-intervals-for-a-proportion","title":"Coverage of interval estimates for a proportion","thumbnail":"eb7eb6ac4017daabc402ae141062f54d36839ef37bd86550c46c2a0f0a4e8662","access_level":"public","create_time":"$D2024-05-28T18:00:37.414Z","update_time":"$D2026-09-03T10:58:31.219Z","view_time":null,"first_public_time":"$D2024-05-28T19:45:10.985Z","fork":{"id":"d8c9f3903311b281","slug":"confidence-interval-for-a-proportion","title":"Confidence interval for a proportion","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":"pro_2024"},"version":971,"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"},{"type":"team","id":"7592a7ab121f74a0","login":"mattiasvillani","name":"Mattias Villani","avatar_url":"https://avatars-next.observableusercontent.com/a440aaf39a3abefb86cbcc758fd7a5ab42555ae68d6825b52ec006f385a94d1d","tier":"pro_2024","role":"viewer"}],"comments":[],"custom_thumbnail":false,"trashed":false,"subscribed":null,"forks":0,"collection":{"id":"3247547e8812c274","slug":"statistics","title":"Basic Statistics","owner":{"id":"7592a7ab121f74a0","login":"mattiasvillani","name":"Mattias Villani"}},"type":"notebook","body":{"title":"Coverage of interval estimates for a proportion","cells":[{"id":0,"value":"# Coverage of interval estimates for a proportion\n","pinned":false,"hidden":false,"mode":"md"},{"id":652,"value":"This notebook investigates the coverage probabilities of some widely used confidence and Bayesian intervals for a proportion ${tex`p`} in the iid Bernoulli model\n${tex.block`X_1,\\ldots,X_N \\mid p \\sim \\mathrm{Bernoulli}(p)`}","pinned":false,"hidden":false,"mode":"md"},{"id":144,"value":"md`\n#### Wald interval\nThe Wald confidence interval with approximate coverage probability ${tex`1-\\alpha`} for a population proportion *p* is\n\n${tex.block`\\hat p \\pm z_{\\alpha/2} \\sqrt{\\frac{\\hat p (1-\\hat p)}{n}}`}\n\nwhere ${tex`\\hat p`} is the proportion in the sample and ${tex`z_{\\alpha/2}`} is obtained from the standard normal distribution as the value that has ${tex`\\alpha/2`} probability to the right of the value. \n\nThe Wald confidence interval is based on a [normal approximation to the binomial distribution](https://observablehq.com/@mattiasvillani/approximating-the-binomial-distribution) which is known to be accurate if both ${tex`np \\geq 10`} and ${tex`n(1-p) \\geq 10`}. \n\n---\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":391,"value":"md`\n#### The Wilson interval\n\nThe Wilson interval uses the formula:\n\n${tex.block`\\tilde p \\pm z_{\\alpha/2} \\frac{\\sqrt{\\frac{\\hat p (1-\\hat p)}{n} + \\frac{z_{\\alpha/2}^2}{4 n^2}}}{1+z_{\\alpha/2}^2/n}`}\n\nwhere the midpoint of the interval is ${tex.block`\\tilde p = \\hat p + \\frac{z_{\\alpha/2}^2}{2n} \\frac{1-2\\hat p}{1+z_{\\alpha/2}^2/n}`}\n\nThe Wilson interval is a little more complicated and less intuitive formula than the Wald interval, but it is often recommended since it tends to have a more accurate coverage^{[[1](#footnote1)]}; its actual coverage is often closer than the nominal coverage, compared to the Wald interval.\n\n---\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":395,"value":"md`\n\n#### Bayesian credible interval\nThe conjugate prior for ${tex`p`} in the iid Bernoulli model is the [Beta distribution](https://observablehq.com/@mattiasvillani/beta-distribution)\n\n${tex.block`p \\sim \\operatorname{Beta}(a, b).`}\n\nCombining this prior with the Bernoulli likelihood gives the posterior\n\n${tex.block`p \\vert x_1,\\ldots,x_n \\sim \\mathrm{Beta}(a+s,a+f)`}\n\nwhere *s* is the number of successes and *f=n-s* is the number of failures.\n\nAn (equal tail) ${tex`1-\\alpha`} [credible interval](https://observablehq.com/@mattiasvillani/bayesian-credible-intervals) is the interval ${tex`[l,u]`} l where *l* and *u* are the 2.5% and 97.5% quantiles in the ${tex`\\mathrm{Beta}(a+s,b+f)`} posterior.\n\nYou can change the prior hyperparameters ${tex`a`} and ${tex`b`} below to see how the coverage of the interval changes.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":412,"value":"viewof priorparam = Inputs.form([\n Inputs.range([0.1, 50], {value: 0.5, label: tex`a`, step: 0.1}),\n Inputs.range([0.1, 50], {value: 0.5, label: tex`b`, step: 0.1}),\n ])","pinned":false,"hidden":false,"mode":"ojs"},{"id":533,"value":"---","pinned":false,"hidden":false,"mode":"md"},{"id":5,"value":"viewof settings = Inputs.form([\n Inputs.range([0, 100], {value: 90, label: \"confidence level\", step: 1}),\n Inputs.range([1, 200], {value: 50, label: \"sample size, n\", step: 1}),\n Inputs.range([0, 1], {value: 0, label: \"lower limit in plot\", step: 0.1}),\n ])","pinned":false,"hidden":false,"mode":"ojs"},{"id":456,"value":"viewof plotselector = Inputs.checkbox([\"wald\", \"wilson\", \"bayes\"], {label: \"plot interval type:\", value: [\"wald\"]})","pinned":false,"hidden":false,"mode":"ojs"},{"id":558,"value":"textinfo = md`\nAverage absolute deviation to target coverage:
\nWald: ${abserrors[0].toFixed(4)} | \nWilson: ${abserrors[1].toFixed(4)} | \nBayes: ${abserrors[2].toFixed(4)} \n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":447,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":808,"value":"viewof plotValidProb = Inputs.toggle({label: md`Plot where normal approx is valid`, value: false})","pinned":false,"hidden":false,"mode":"ojs"},{"id":894,"value":"

\n#### Footnotes","pinned":false,"hidden":false,"mode":"md"},{"id":887,"value":"footnote1 = md`1. [^](#paragraph1) \nAndersson, P. G. (2023). The Wald confidence interval for a binomial p as an illuminating “bad” example, *The American Statistician*, 77:4, 443-448.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":12,"value":"
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n\n---\n## 🙈 Appendix 🙈\n\nHere is just code and details that you probably do not want to read.","pinned":false,"hidden":false,"mode":"md"},{"id":246,"value":"#### Imports","pinned":false,"hidden":false,"mode":"md"},{"id":15,"value":"jstat = require('jstat')","pinned":false,"hidden":false,"mode":"ojs"},{"id":206,"value":"import {Scrubber} from \"@mbostock/scrubber\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":248,"value":"#### Functions","pinned":false,"hidden":false,"mode":"md"},{"id":329,"value":"wald = function(s, n, zcrit){\n const p̂ = s/n\n const se = Math.sqrt( p̂*(1 - p̂)/n)\n return [p̂ - zcrit*se, p̂ + zcrit*se]\n} ","pinned":false,"hidden":false,"mode":"ojs"},{"id":339,"value":"wilson = function(s, n, zcrit){\n const p̂ = s/n\n const ptilde = p̂ + (zcrit**2/(2*n))*((1-2*p̂)/(1+zcrit**2/n))\n const se = Math.sqrt( p̂*(1 - p̂)/n + (zcrit**2/(4*n**2)) )/(1+zcrit**2/n) \n return [ptilde - zcrit*se, ptilde + zcrit*se]\n} ","pinned":false,"hidden":false,"mode":"ojs"},{"id":422,"value":"bayesian = function(s, n, confidence, alpha, beta){\n const lower = jstat.beta.inv((1 - confidence/100)/2, alpha + s, beta + n - s)\n const upper = jstat.beta.inv(1 - (1 - confidence/100)/2, alpha + s, beta + n - s)\n return [lower, upper]\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":250,"value":"#### Data","pinned":false,"hidden":false,"mode":"md"},{"id":462,"value":"colors = [\"#6C8EBF\", \"#c0a34d\", \"#780000\", \"#007878\", \"#bf9d6c\", \"#3A6B35\", \"#b5c6df\",\"#eadaaa\", \"#bb989a\", \"#98bbb9\", \"#bf8d6c\", \"#CBD18F\"]","pinned":false,"hidden":false,"mode":"ojs"},{"id":28,"value":"confidence = settings[0]","pinned":false,"hidden":false,"mode":"ojs"},{"id":26,"value":"n = settings[1]","pinned":false,"hidden":false,"mode":"ojs"},{"id":488,"value":"lower_y = settings[2]","pinned":false,"hidden":false,"mode":"ojs"},{"id":336,"value":"zcrit = jstat.normal.inv(1 - (1 - confidence/100)/2, 0 , 1)","pinned":false,"hidden":false,"mode":"ojs"},{"id":417,"value":"alpha = priorparam[0]","pinned":false,"hidden":false,"mode":"ojs"},{"id":419,"value":"beta = priorparam[1]","pinned":false,"hidden":false,"mode":"ojs"},{"id":683,"value":"lower = 0.001// onlyValid_p ? 10/n : 0.001","pinned":false,"hidden":false,"mode":"ojs"},{"id":686,"value":"upper = 0.999 // onlyValid_p ? 1-10/n : 0.999","pinned":false,"hidden":false,"mode":"ojs"},{"id":688,"value":"pgrid = d3.range(lower, upper, (upper-lower)/1000)","pinned":false,"hidden":false,"mode":"ojs"},{"id":438,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":550,"value":"abserrors = [d3.mean(coveragedata.map(d => Math.abs(d.wald - confidence/100))), \n d3.mean(coveragedata.map(d => Math.abs(d.wilson - confidence/100))), \n d3.mean(coveragedata.map(d => Math.abs(d.bayes - confidence/100)))\n]","pinned":false,"hidden":false,"mode":"ojs"},{"id":671,"value":"priorcoverage = d3.sum(coveragedata.map(d => d.priorpdf*d.bayes*(pgrid[1]-pgrid[0])))","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W34","stdlib":"1","titling":"auto"}} 30:"$2f:body"