31:T4c1,{ const svg = d3.select(DOM.svg(size.width, size.height)); svg.append("rect") .attr("fill", "#f0f0f0") .attr("x", margin.left) .attr("y", margin.top) .attr("width", size.width - margin.left - margin.right) .attr("height", size.height - margin.top - margin.bottom); svg.append('g').call(xGrid); svg.append('g').call(yGrid); svg.append('g').call(xAxis); svg.append("g").call(yAxis); svg.append("path") .attr("fill", "#ccc") .attr("fill-opacity", 0.5) .attr("stroke", "none") .attr("d", area(intervalData)); svg.append("line") .datum([ [x.domain()[0], predict(x.domain()[0])], [x.domain()[1], predict(x.domain()[1])] ]) .attr("stroke", "#333") .attr("stroke-width", 1.5) .attr("x1", d => x(d[0][0])) .attr("x2", d => x(d[1][0])) .attr("y1", d => y(d[0][1])) .attr("y2", d => y(d[1][1])); svg.append("g") .attr("stroke", "white") .attr("stroke-width", 1.5) .attr("fill", "steelblue") .selectAll("circle") .data(data) .join("circle") .attr("cx", d => x(d.x)) .attr("cy", d => y(d.y)) .attr("r", 4) .append("title").text(d => d.name); return svg.node(); }32:T566,md` The chart above shows the relationship between engine power (*y*-axis) and fuel consumption (*x*-axis) for the ${data.length} cars of the *mtcars* dataset, along with a [linear regression](https://en.wikipedia.org/wiki/Linear_regression) line. The grey bands around the regression line represent the range in which the *true* regression line lies at a certain level of confidence. They visualize all [confidence intervals](https://en.wikipedia.org/wiki/Confidence_interval) for every possible ${tex`x`} and are tightest where the data is grouped more densely. The confidence interval around the regression line for a given ${tex`x`} is calculated using: ${tex.block` ŷ_{ci} = ŷ \pm t_{n-2,\alpha/2} \sqrt{\frac{\sum_{i=1}^n (y_i - ŷ)²}{n-2}} \sqrt{\frac{1}{n} + \frac{(x - \bar{x})²}{\sum_{i=1}^{n}(x_i - \bar{x})²}} `} where the critical [${tex`t`}-value](https://en.wikipedia.org/wiki/Student%27s_t-test) is ${tex`t_{n-2,\alpha/2}`} with ${tex`n`} − 2 degrees of freedom. It is calculated using Hill's approximated inverse t-distribution in this example.

Hill, G. W. (1970), Algorithm 396: Student's t-quantiles. *Communications of the ACM* 13(10), 619–620.
Wichura, M. J. (1988), Algorithm AS 241: The percentage points of the normal distribution. *Applied Statistics* 37, 477–484.

`33:T61e,/* Hill's approximated inverse t-distribution calculates t given df and two-tail probability: Hill, G. W. (1970), Algorithm 396: Student's t-quantiles. Communications of the ACM, 13(10), 619–620. Result should be "correct to at least 6 significant digits even for the analytic continuation through noninteger values of n > 5". For higher precision (used in R, ...) see: Hill, G. W. (1981) Remark on Algorithm 396, ACM Transactions on Mathematical Software, 7, 250–1. Excel: TINV(); R: qt() */ function inverseT(p, df) { const { sin, cos, sqrt, pow, exp, PI} = Math; let a, b, c, d, t, x, y; if (df == 1) return cos(p * PI / 2) / sin( p * PI / 2); if (df == 2) return sqrt(2 / (p * (2 - p)) - 2); a = 1 / (df - 0.5); b = 48 / (a * a); c = ((20700 * a / b - 98) * a - 16) * a + 96.36; d = ((94.5 / (b + c) - 3) / b + 1) * sqrt(a * PI * 0.5) * df; x = d * p; y = pow(x, 2 / df); if (y > 0.05 + a) { // The procedure normdev(p) is assumed to return a negative normal // deviate at the lower tail probability level p, e.g. -2.32 for p = 0.01. x = normdev(p / 2); y = x * x; if (df < 5) c = c + 0.3 * (df - 4.5) * (x + 0.6); c = (((0.05 * d * x - 5) * x - 7) * x - 2) * x + b + c; y = (((((0.4 * y + 6.3) * y + 36) * y + 94.5) / c - y - 3) / b + 1) * x; y = a * y * y; if (y > 0.002) y = exp(y) - 1; else y = 0.5 * y * y + y; } else { y = ((1 / (((df + 6) / (df * y) - 0.089 * d - 0.822) * (df + 2) * 3) + 0.5 / (df + 4)) * y - 1) * (df + 1) / (df + 2) + 1 / y; } return sqrt(df * y); }34:Tc62,/** Returns the normal deviate Z corresponding to a given lower tail area of P; Z is accurate to about 1 part in 10**16. Wichura, M. J. (1988) Algorithm AS 241: The percentage points of the normal distribution. Applied Statistics, 37, 477–484. extension of: Beasley, J. D./ Springer, S. G. (1977), The percentage points of the NormalDistribution, Applied Statistics. 26, 118–121. Excel: STANDNORMINV(); R: qnorm() **/ function normdev(p) { if (p < 0 || p > 1) return false; if (p == 0) return -Infinity; if (p == 1) return Infinity; let a0 = 3.3871328727963666080e+0, a1 = 1.3314166789178437745e+2, a2 = 1.9715909503065514427e+3, a3 = 1.3731693765509461125e+4, a4 = 4.5921953931549871457e+4, a5 = 6.7265770927008700853e+4, a6 = 3.3430575583588128105e+4, a7 = 2.5090809287301226727e+3, b1 = 4.2313330701600911252e+1, b2 = 6.8718700749205790830e+2, b3 = 5.3941960214247511077e+3, b4 = 2.1213794301586595867e+4, b5 = 3.9307895800092710610e+4, b6 = 2.8729085735721942674e+4, b7 = 5.2264952788528545610e+3, c0 = 1.42343711074968357734e+0, c1 = 4.63033784615654529590e+0, c2 = 5.76949722146069140550e+0, c3 = 3.64784832476320460504e+0, c4 = 1.27045825245236838258e+0, c5 = 2.41780725177450611770e-1, c6 = 2.27238449892691845833e-2, c7 = 7.74545014278341407640e-4, d1 = 2.05319162663775882187e+0, d2 = 1.67638483018380384940e+0, d3 = 6.89767334985100004550e-1, d4 = 1.48103976427480074590e-1, d5 = 1.51986665636164571966e-2, d6 = 5.47593808499534494600e-4, d7 = 1.05075007164441684324e-9, e0 = 6.65790464350110377720e+0, e1 = 5.46378491116411436990e+0, e2 = 1.78482653991729133580e+0, e3 = 2.96560571828504891230e-1, e4 = 2.65321895265761230930e-2, e5 = 1.24266094738807843860e-3, e6 = 2.71155556874348757815e-5, e7 = 2.01033439929228813265e-7, f1 = 5.99832206555887937690e-1, f2 = 1.36929880922735805310e-1, f3 = 1.48753612908506148525e-2, f4 = 7.86869131145613259100e-4, f5 = 1.84631831751005468180e-5, f6 = 1.42151175831644588870e-7, f7 = 2.04426310338993978564e-15; let q = p - 0.5, r, z; // p close to 0.5 if (Math.abs(q) <= 0.425) { r = 0.180625 - q * q; z = q * (((((((a7 * r + a6) * r + a5) * r + a4) * r + a3) * r + a2) * r + a1) * r + a0) / (((((((b7 * r + b6) * r + b5) * r + b4) * r + b3) * r + b2) * r + b1) * r + 1); return z; } if (q > 0) r = 1 - p; else r = p; r = Math.sqrt(-Math.log(r)); // p neither close to 0.5 nor 0 or 1 if (r <= 5) { r += -1.6; z = (((((((c7 * r + c6) * r + c5) * r + c4) * r + c3) * r + c2) * r + c1) * r + c0) / (((((((d7 * r + d6) * r + d5) * r + d4) * r + d3) * r + d2) * r + d1) * r + 1); } // p near 0 or 1 else { r += -5; z = (((((((e7 * r + e6) * r + e5) * r + e4) * r + e3) * r + e2) * r + e1) * r + e0) / (((((((f7 * r + f6) * r + f5) * r + f4) * r + f3) * r + f2) * r + f1) * r + 1); } if (q < 0.0) z = -z; return z; }35:[] 2f:{"title":"Linear regression with confidence bands","cells":[{"id":0,"value":"md`# Linear regression with confidence bands`","pinned":false,"hidden":false,"mode":"ojs"},{"id":98,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":630,"value":"viewof confidenceLevel = radio({\n title: 'Confidence level',\n options: [\n { label: \"95%\", value: 0.05 },\n { label: \"96%\", value: 0.04 },\n { label: \"98%\", value: 0.02 },\n { label: \"99%\", value: 0.01 },\n { label: \"99.99%\", value: 0.0001 },\n ],\n value: 0.05\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":682,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":451,"value":"md`---\n## Appendix`","pinned":false,"hidden":false,"mode":"ojs"},{"id":190,"value":"predict = linearRegression(data)","pinned":false,"hidden":false,"mode":"ojs"},{"id":16,"value":"interval = confidenceInterval(data, predict)","pinned":false,"hidden":false,"mode":"ojs"},{"id":578,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":893,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":4,"value":"function linearRegression(data) {\n let n = data.length, sumX = 0, sumY = 0, sumXY = 0, sumX2 = 0;\n for (let i = 0; i < data.length; ++i) {\n const xi = data[i].x;\n const yi = data[i].y;\n if (!Number.isFinite(xi) || !Number.isFinite(yi)) { --n; continue; }\n sumX += xi;\n sumY += yi;\n sumXY += xi * yi;\n sumX2 += xi * xi;\n }\n const slope = (n * sumXY - sumX * sumY) / (n * sumX2 - sumX * sumX);\n const intercept = (sumY - slope * sumX) / n;\n return x => slope * x + intercept;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":42,"value":"// https://stats.stackexchange.com/questions/101318/understanding-shape-and-calculation-of-confidence-bands-in-linear-regression\nfunction confidenceInterval(data, predict) {\n const mean = d3.sum(data, d => d.x) / data.length;\n let a = 0, b = 0;\n for (let i = 0; i < data.length; ++i) {\n a += Math.pow(data[i].x - mean, 2);\n b += Math.pow(data[i].y - predict(data[i].x), 2);\n }\n const sy = Math.sqrt(b / (data.length - 2));\n const t = inverseT(+confidenceLevel, data.length - 2) \n return function(x) {\n const Y = predict(x);\n const se = sy * Math.sqrt(1 / data.length + Math.pow(x - mean, 2) / a);\n return { x, left: Y - t * se, right: Y + t * se };\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":145,"value":"intervalData = {\n const domain = x.domain();\n const step = (domain[1] - domain[0]) / 100;\n return d3.range(domain[0], domain[1] + step, step).map(d => interval(d));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":380,"value":"data = mtcars.map(car => ({ name: car.name, x: 1 / car.mpg, y: car.hp}))","pinned":false,"hidden":false,"mode":"ojs"},{"id":454,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs"},{"id":332,"value":"area = d3.area()\n .x(d => x(d.x))\n .y0(d => y(d.left))\n .y1(d => y(d.right))","pinned":false,"hidden":false,"mode":"ojs"},{"id":106,"value":"x = d3.scaleLinear()\n .domain(d3.extent(data, d => d.x))\n\t.range([margin.left, size.width - margin.right])\n .nice()","pinned":false,"hidden":false,"mode":"ojs"},{"id":108,"value":"y = d3.scaleLinear()\n\t.domain([d3.min(intervalData, d => d.left), d3.max(intervalData, d => d.right)])\n\t.range([size.height - margin.bottom, margin.top])\n .nice()","pinned":false,"hidden":false,"mode":"ojs"},{"id":803,"value":"xGrid = g => g\n .attr(\"transform\", `translate(0,${size.height - margin.bottom})`)\n .call(d3.axisBottom(x).tickSize(-size.height + margin.top + margin.bottom))\n .call(g => g.select(\".domain\").remove())\n .call(g => g.selectAll(\"text\").remove())\n .call(g => g.selectAll(\"line\").attr(\"stroke\", \"white\"));","pinned":false,"hidden":false,"mode":"ojs"},{"id":814,"value":"yGrid = g => g\n .attr(\"transform\", `translate(${margin.left},0)`)\n .call(d3.axisLeft(y).tickSize(-size.width + margin.left + margin.right))\n .call(g => g.select(\".domain\").remove())\n .call(g => g.selectAll(\"text\").remove())\n .call(g => g.selectAll(\"line\").attr(\"stroke\", \"white\"));","pinned":false,"hidden":false,"mode":"ojs"},{"id":111,"value":"xAxis = g => g\n .attr(\"transform\", `translate(0,${size.height - margin.bottom})`)\n .call(d3.axisBottom(x))\n .call(g => g.select(\".domain\").remove())\n .call(g => g.append(\"text\")\n .attr(\"x\", size.width - margin.right - 4)\n .attr(\"y\", -4)\n .attr(\"fill\", \"#000\")\n .attr(\"font-weight\", \"bold\")\n .attr(\"text-anchor\", \"end\")\n .text(\"Gallons p. mile\"))","pinned":false,"hidden":false,"mode":"ojs"},{"id":113,"value":"yAxis = g => g\n .attr(\"transform\", `translate(${margin.left},0)`)\n .call(d3.axisLeft(y))\n .call(g => g.select(\".domain\").remove())\n .call(g => g.select(\".tick:last-of-type text\").clone()\n .attr(\"x\", 4)\n .attr(\"y\", 8)\n .attr(\"text-anchor\", \"start\")\n .attr(\"font-weight\", \"bold\")\n .text(\"Horsepower\"))","pinned":false,"hidden":false,"mode":"ojs"},{"id":103,"value":"margin = ({top: 12, right: 20, bottom: 26, left: 26})","pinned":false,"hidden":false,"mode":"ojs"},{"id":99,"value":"size = ({width, height: width / 3 * 2})","pinned":false,"hidden":false,"mode":"ojs"},{"id":458,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs"},{"id":367,"value":"mtcars = require(\"@observablehq/mtcars\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":76,"value":"d3 = require(\"d3@5\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":632,"value":"import { radio } from \"@jashkenas/inputs\"","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W35","stdlib":"1","titling":"auto"}