2f:T552,# Curve fitting with B-splines B-splines are a powerful tool commonly used in statistics to model smooth functions. They can be used in settings where we wish to fit a curve to data without making strong functional form assumptions. We will call these curves _B-spline functions._ B-spline functions are defined in terms of a set of basis functions, so we'll start by defining what a B-spline basis function is. ## B-spline Basis Functions We need two ingredients to build a B-spline basis function: the order ${tex`d`} and a sequence of non-decreasing knots ${tex`t_1, t_2, \dots, t_k`}. As a running example, let's use ${tex`[${knots}]`} as a series of knots (so ${tex`k=5`} in this example.) ${sidenote(md`Note: you also see B-splines referred to by their _degree_. The degree ${tex`p`} of a B-spline is defined as its order minus one: ${tex`p = d - 1`}.`)} Let ${tex`B_{i,d}`} be the ${tex`i^{th}`} B-spline of order ${tex`d`}. These functions are defined recursively, starting with ${tex`d=1`}, in which case the function is piecewise constant: ${tex.block` B_{i,1}(x) = \begin{cases} 1, & t_i \leq x < t_{i+1},\\ 0, & \mathrm{otherwise.} \end{cases} `} The plot below shows the 1st order basis functions for our set of example knots. There's not much to see here yet, as the basis functions are just piecewise constant between each knot location.30:T72e,One little thing to notice here is that at ${tex`x=1`}, the last basis function evaluates to ${tex`0`} (this is why you see the vertical line right at ${tex`x =1`} on the plot.) This is a consequence of how we defined the basis functions to be equal to ${tex`1`} only if ${tex`t_i \leq x < t_{i+1}`}. The strictly less than sign means that the basis functions will always be ${tex`0`} at the last knot location. We could work around this by redefining the final basis function to be ${tex`1`} at the last knot location, but for now for simplicity's sake we'll leave it as is. Now that we have defined a base case, we can define higher order B-spline functions via a recursive formula: ${tex.block` B_{i,d}(x) = \omega_{i,d} B_{i,d-1}(x) + (1 - \omega_{i,d}) B_{i + 1, d - 1}(x) `} where ${tex.block` \omega_{i,d}(x) = \begin{cases} \frac{x - t_i}{t_{i+d-1} - t_i}, & t_i \neq t_{i + d - 1} \\ 0, & \mathrm{otherwise}. \end{cases} `} There is one little hiccup in this formula: notice how it references a knot at ${tex`t_{i+d-1}`}. This overruns the number of knots that we've defined so far. For this formula to be well defined we need to add ${tex`d-1`} new knots to the start and end of our knot vector, creating a set of extended knots: ${tex.block` \underset{d-1}{\underbrace{t_1, t_1, t_1}}, \underset{\mathrm{original\ knots}}{\underbrace{t_1, t_2, \dots, t_k}}, \underset{d-1}{\underbrace{t_k, t_k, t_k}} `} We also need to define more basis functions than before in order to cover all of the space between the first and last knot. Specifically, we need a total of ${tex`k + d - 2`} basis functions, where ${tex`k`} is the original number of knots and ${tex`d`} is the spline order. The plot below shows the basis functions for different spline orders: you can move the slider to change the order and the plot will update. 31:T50a,A couple interesting things to note from these plots: - Each basis function is only positive over a subset of the domain. As the order increases, the basis functions are positive over longer intervals. - As you increase the order of the basis functions, they overlap one another more. These two things help explain how B-spline functions, which combine the basis functions together, become more smooth as the B-spline order increases. We'll talk about how B-spline functions work next. ## B-spline Functions So far we've described the definition of B-spline basis functions. Next we will look at how to use this as a foundation to define a B-spline function. As before, suppose we have knots ${tex`t_1, t_2, \dots, t_k`} (and a set of extended knots, constructed as before.) A B-spline function ${tex`f(x)`} of order ${tex`d`} is defined as the linear combination of ${tex`k+d-2`} basis functions: ${tex.block` f(x) = \sum_{i=1}^{k + d - 2} \alpha_{i} B_{i,d}(x). `} The parameters ${tex`\alpha_1, \alpha_1, \alpha_{k+d-2}`} are called _spline coefficients_. The plot below shows the spline basis functions, like before, but now you can click and drag the circles to change the coefficients associated with each basis. The second plot shows the resulting spline function ${tex`f(x)`}.32:T970,plot2 = { const height = plot_height; const plot_width = width / 2 - 10; const svg = d3.select(DOM.svg(plot_width, height)) .attr("viewBox", [0, 0, plot_width, height]); const margin = { top: 20, right: 50, bottom: 50, left: 100 }; const node = svg.node(); const ylim = [0, 3]; const B_scaled = B.map((row, i) => row.map((d) => d * spline_coefs[i])); const knot_data = []; mutable spline_coefs; for(let i = 0; i < spline_coefs.length; i++) { const spline_max_index = indexOfMax(B[i]); indexOfMax(B_scaled[i].map((d) => -d)); knot_data[i] = ({ x: grid[spline_max_index], y: spline_coefs[i], index: i }); } const x = d3.scaleLinear() .domain([0, 1]) .range([margin.left, plot_width - margin.right]) const y = d3.scaleLinear() .domain(ylim) .range([height - margin.bottom, margin.top]); svg.append("g") .attr("class", "axis axis--x") .attr("transform", `translate(${-0.5},${height - margin.bottom})`) .call(d3.axisBottom(x)); svg.append("g") .attr("class", "axis axis--y") .attr("transform", `translate(${margin.left - 30},0)`) .call(d3.axisLeft(y).ticks(5)); svg.append("text") .style("font-style", "italic") .text("x") .attr("transform", `translate(${plot_width/2}, ${height})`); svg.append("text") .style("font-style", "italic") .attr("text-anchor", "end") .text("f(x)") .attr("transform", `translate(${margin.left - 70}, ${margin.top + 10})`); const line = d3.line() .x((d) => x(d.x)) .y((d) => y(d.y)); svg.selectAll('path.f') .data([spline_ys]) .enter() .append('path') .attr('class', 'f') .attr('d', (d) => line(d)) .attr('stroke', "black") .attr('fill', 'transparent') const drag = d3.drag() .container(svg) .clickDistance(100) .on("drag", function(event, d) { spline_coefs[d.index] = Math.min(ylim[1], Math.max(ylim[0], y.invert(event.sourceEvent.offsetY))); mutable spline_coefs = Array.from(spline_coefs); }); const accent = d3.scaleOrdinal(d3.schemeCategory10); svg.selectAll('circle.knot') .data(knot_data) .enter() .append('circle') .attr('cx', (d) => x(d.x)) .attr('cy', (d) => y(d.y)) .attr('r', 6) .attr('fill', (d, i) => accent(i)) .attr('cursor', 'move') .call(drag); return node; }33:T4e5,## Odds and ends - A _cardinal B-spline_ is a B-spline with knots that are placed a constant width apart. All the splines shown in this notebook are cardinal splines, because the knots are spaced evenly between ${tex`0`} and ${tex`1`}. - If you want to write a computer program to evaluate a spline function at some arbitrary ${tex`x`}, you can do better than calculating ${tex` f(x) = \sum_{i=0}^{k + d} \alpha_{i} B_{i,d}(x) `} directly. At any one ${tex`x`}, most of the spline basis functions will evaluate to zero, and won't contribute anything to the result, so it's a waste of time calculating them. A better option is [De Boor's Algorithm](https://en.wikipedia.org/wiki/De_Boor%27s_algorithm), which is smart enough to only calculate the basis functions which contribute to the result. - A spline function of order ${tex`d`} has ${tex`d - 1`} well defined derivatives. In fact, the derivative of a spline function of order ${tex`d`} is itself a spline function of order ${tex`d - 1`}. ## Works Cited I drew heavily from the [Splines in Stan](https://mc-stan.org/users/documentation/case-studies/splines_in_stan.html) case study by Milad Kharratzadeh. If you're interested in using B-splines in Stan, that case study is a great place to start.34:Tee2,plot_basis_functions_interactive = (grid, B, order, knots, coefs, x_title) => { const height = plot_height; const plot_width = width / 2 - 10; const svg = d3.select(DOM.svg(plot_width, height)) .attr("viewBox", [0, 0, plot_width, height]); const margin = { top: 20, right: 50, bottom: 50, left: 60 }; const node = svg.node(); const ylim = [0, 3]; const B_scaled = B.map((row, i) => row.map((d) => d * spline_coefs[i])); const knot_data = []; mutable spline_coefs; for(let i = 0; i < spline_coefs.length; i++) { const spline_max_index = indexOfMax(B[i]); knot_data[i] = ({ x: grid[spline_max_index], y: spline_coefs[i], //B_scaled[i][spline_max_index], index: i }); } const x = d3.scaleLinear() .domain([0, 1]) .range([margin.left, plot_width - margin.right]) const y = d3.scaleLinear() .domain(ylim) .range([height - margin.bottom, margin.top]); svg.append("g") .attr("class", "axis axis--x") .attr("transform", `translate(${-0.5},${height - margin.bottom})`) .call(d3.axisBottom(x)); svg.append("g") .attr("class", "axis axis--y") .attr("transform", `translate(${margin.left - 30},0)`) .call(d3.axisLeft(y).ticks(5)); svg.append("text") .style("font-style", "italic") .text("x") .attr("transform", `translate(${width/2}, ${height})`); svg.append("foreignObject") .style("font-style", "italic") .attr("text-anchor", "end") .attr("width", 200) .attr("height", 100) .html(x_title.outerHTML) .attr("transform", `translate(${margin.left - 100}, ${margin.top})`); // Knot location labels const knotLabels = svg.selectAll("g.knot") .data(knots) .enter() .append('g') .attr('class', 'knot') .attr('transform', (d) => `translate(${-0.5 + x(d)}, ${height - margin.bottom + 15})`); knotLabels.append("line") .attr("x1", 0) .attr("y1", 0) .attr("x2", 0) .attr("y2", -15) .attr("stroke-dasharray", "2 2") .attr("stroke", "gray"); /* knotLabels.append("foreignObject") .attr('width', 20) .attr('height', 20) .attr("transform", "translate(-5, 0)") .attr("color", "gray") .html((d, i) => tex`t_{${i + 1}}`.outerHTML) */ const line = d3.line() .x((d, i) => x(grid[i])) .y((d) => y(d)); const accent = d3.scaleOrdinal(d3.schemeCategory10); svg.selectAll('path.basis') .data(B_scaled) .enter() .append('path') .attr('class', 'basis') .attr('d', (d) => line(d)) .attr('stroke', (d, i) => accent(i)) .attr('stroke-dasharray', '2 2') .attr('fill', 'transparent'); /*svg.selectAll('foreignObject.label') .data(B_scaled) .enter() .append('foreignObject') .attr('class', 'label') .attr('width', 115) .attr('height', 30) .attr('transform', (d, i) => { const spliney = spline_coefs[i]; const splinex = grid[indexOfMax(B[i])]; const offset = spline_coefs[i] > 0 ? -30 : 5; return `translate(${x(splinex) - 20}, ${y(spliney) + offset}) scale(0.8)` }) .html((d, i) => tex`\alpha_{${i + 1}} = ${round(coefs[i])}`.outerHTML) .attr('color', (d, i) => accent(i)) .style('pointer-events', 'none');*/ const drag = d3.drag() .container(svg) .clickDistance(100) .on("drag", function(event, d) { spline_coefs[d.index] = Math.min(ylim[1], Math.max(ylim[0], y.invert(event.sourceEvent.offsetY))); mutable spline_coefs = Array.from(spline_coefs); }); svg.selectAll('circle.knot') .data(knot_data) .enter() .append('circle') .attr('cx', (d) => x(d.x)) .attr('cy', (d) => y(d.y)) .attr('r', 6) .attr('fill', (d, i) => accent(i)) .style('cursor', 'move') .call(drag); return node; }35:Ta4e,plot_basis_functions = (grid, B, order, knots, x_title) => { const height = 250; const plot_width = width / 2 - 10; const svg = d3.select(DOM.svg(plot_width, height)) .attr("viewBox", [0, 0, plot_width, height]); const margin = { top: 20, right: 50, bottom: 50, left: 120 }; const node = svg.node(); const ylim = [0, 1.1]; const x = d3.scaleLinear() .domain([0, 1]) .range([margin.left, plot_width - margin.right]) const y = d3.scaleLinear() .domain(ylim) .range([height - margin.bottom, margin.top]); svg.append("g") .attr("class", "axis axis--x") .attr("transform", `translate(${-0.5},${height - margin.bottom})`) .call(d3.axisBottom(x)); svg.append("g") .attr("class", "axis axis--y") .attr("transform", `translate(${margin.left - 10},0)`) .call(d3.axisLeft(y)); svg.append("text") .style("font-style", "italic") .text("x") .attr("transform", `translate(${plot_width/2}, ${height})`); svg.append("foreignObject") .style("font-style", "italic") .attr("text-anchor", "end") .attr("width", 200) .attr("height", 100) .html(x_title.outerHTML) .attr("transform", `translate(${margin.left - 100}, ${margin.top})`); // Knot location labels const knotLabels = svg.selectAll("g.knot") .data(knots) .enter() .append('g') .attr('class', 'knot') .attr('transform', (d) => `translate(${-0.5 + x(d)}, ${height - margin.bottom + 15})`); knotLabels.append("line") .attr("x1", 0) .attr("y1", 0) .attr("x2", 0) .attr("y2", -15) .attr("stroke-dasharray", "2 2") .attr("stroke", "gray"); knotLabels.append("foreignObject") .attr('width', 20) .attr('height', 25) .attr("transform", "translate(-5, 0)") .attr("color", "gray") .html((d, i) => tex`t_{${i + 1}}`.outerHTML) const line = d3.line() .x((d, i) => x(grid[i])) .y((d) => y(d)); let accent = d3.scaleOrdinal(d3.schemeCategory10); svg.selectAll('path.basis') .data(B) .enter() .append('path') .attr('class', 'basis') .attr('d', (d) => line(d)) .attr('stroke', (d, i) => accent(i)) .attr('stroke-dasharray', '2 2') .attr('fill', 'transparent'); svg.selectAll('foreignObject.label') .data(B) .enter() .append('foreignObject') .attr('class', 'label') .attr('width', 35) .attr('height', 30) .attr('transform', (d) => `translate(${x(grid[indexOfMax(d)])}, ${y(d3.max(d)) - 30})`) .html((d, i) => tex`B_{${i + 1}, ${order}}(x)`.outerHTML) .attr('color', (d, i) => accent(i)); return node; }36:T5ab,{ return null; const height = 250; const svg = d3.select(DOM.svg(width, height)) .attr("viewBox", [0, 0, width, height]); const margin = { top: 20, right: 50, bottom: 50, left: 120 }; const node = svg.node(); const ylim = [-5, 5]; //d3.extent(adjusted_B.flat()); const x = d3.scaleLinear() .domain([0, 1]) .range([margin.left, width - margin.right]) const y = d3.scaleLinear() .domain(ylim) .range([height - margin.bottom, margin.top]); svg.append("g") .attr("class", "axis axis--x") .attr("transform", `translate(${-0.5},${height - margin.bottom})`) .call(d3.axisBottom(x)); svg.append("g") .attr("class", "axis axis--y") .attr("transform", `translate(${margin.left - 10},0)`) .call(d3.axisLeft(y)); svg.append("text") .style("font-style", "italic") .text("x") .attr("transform", `translate(${width/2}, ${height})`); svg.append("text") .style("font-style", "italic") .attr("text-anchor", "end") .text("B-spline basis") .attr("transform", `translate(${margin.left - 50}, ${margin.top + 10}) rotate(-90)`); const line = d3.line() .x((d, i) => x(grid[i])) .y((d) => y(d)); svg.selectAll('path.basis') .data(adjusted_B) .enter() .append('path') .attr('class', 'basis') .attr('d', (d) => line(d)) .attr('stroke', "black") .attr('fill', 'transparent'); return node; }37:[] 2d:{"title":"Curve fitting with B-splines","cells":[{"id":0,"value":"$2f","pinned":false,"hidden":false,"mode":"md"},{"id":222,"value":"plot_basis_functions(grid, B_order0, 1, knots, tex`B_{i,1}(x)`)","pinned":false,"hidden":false,"mode":"ojs"},{"id":199,"value":"$30","pinned":false,"hidden":false,"mode":"md"},{"id":305,"value":"Inputs.bind(Inputs.range([1, 5], { value: 4, step: 1, label: \"Spline order\" }), viewof spline_degree)","pinned":false,"hidden":false,"mode":"ojs"},{"id":303,"value":"plot_basis_functions(grid, B, spline_order, knots, tex`B_{i,${spline_order}}(x)`)","pinned":false,"hidden":false,"mode":"ojs"},{"id":329,"value":"$31","pinned":false,"hidden":false,"mode":"md"},{"id":10,"value":"viewof num_knots = Inputs.range([2, 8], { value: 5, step: 1, label: \"Number of knots\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":85,"value":"spline_degree_view = Inputs.bind(Inputs.range([1, 6], { value: 3, step: 1, label: \"Spline degree\" }), viewof spline_degree)","pinned":false,"hidden":false,"mode":"ojs"},{"id":692,"value":"plots = html`\n${plot1} ${plot2}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":218,"value":"plot1 = plot_basis_functions_interactive(grid, B, spline_order, knots, spline_coefs, \"B-spline basis\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":4,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":720,"value":"params = new URL(document.location).searchParams","pinned":true,"hidden":false,"mode":"ojs"},{"id":718,"value":"plot_height = params.get(\"height\") == null ? 250 : parseInt(params.get(\"height\") - 100)","pinned":true,"hidden":false,"mode":"ojs"},{"id":423,"value":"$33","pinned":false,"hidden":false,"mode":"md"},{"id":625,"value":"sidenote = (content) => html`

\n${content}\n

`","pinned":false,"hidden":false,"mode":"ojs"},{"id":474,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":133,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":481,"value":"round = d3.format(\"0.2r\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":161,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":26,"value":"knots = {\n const k = [...Array(num_knots)].map((x, i) => i / (num_knots - 1));\n return k;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":87,"value":"num_basis = num_knots + spline_order - 2","pinned":false,"hidden":false,"mode":"ojs"},{"id":52,"value":"build_ext_knots = (knots, degree) => [Array(degree).fill(knots[0]), knots, Array(degree).fill(knots[knots.length - 1])].flat()","pinned":false,"hidden":false,"mode":"ojs"},{"id":23,"value":"mutable spline_coefs = Array(num_basis).fill(1)","pinned":false,"hidden":false,"mode":"ojs"},{"id":61,"value":"grid = [...Array(101)].map((x, i) => i / 100).concat(knots.filter((d) => d != 0).map((d) => d - 0.0001)).concat(knots.map((d) => d + 0.0001)).sort()","pinned":false,"hidden":false,"mode":"ojs"},{"id":63,"value":"spline_ys = spline_f(grid, knots, spline_coefs, spline_order - 1)","pinned":false,"hidden":false,"mode":"ojs"},{"id":45,"value":"spline_f = (x, knots, coefs, degree) => {\n const B = basis_functions(x, knots, degree);\n\n const y = mathjs.multiply(coefs, B);\n //const y = 0;\n return x.map((x, i) => ({\n x: x,\n y: y._data[i]\n }));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":128,"value":"B = basis_functions(grid, knots, spline_order - 1)._data","pinned":false,"hidden":false,"mode":"ojs"},{"id":220,"value":"B_order0 = basis_functions(grid, knots, 0)._data","pinned":false,"hidden":false,"mode":"ojs"},{"id":146,"value":"adjusted_B = B.map((row, index) => row.map((x) => x * spline_coefs[index]))","pinned":false,"hidden":false,"mode":"ojs"},{"id":566,"value":"viewof spline_degree = Inputs.input(2)","pinned":true,"hidden":false,"mode":"ojs"},{"id":684,"value":"spline_order = spline_degree + 1","pinned":true,"hidden":false,"mode":"ojs"},{"id":124,"value":"basis_functions = (x, knots, degree) => {\n const order = degree + 1;\n const num_basis = knots.length + degree - 1;\n const ext_knots = build_ext_knots(knots, degree);\n const mat = [];\n for(let index = 0; index < num_basis; index++) {\n mat[index] = build_b_spline(x, ext_knots, index, order);\n }\n //mat[mat.length] = Array(x.length).fill(0);\n\n const B = mathjs.matrix(mat);\n\n return B;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":42,"value":"build_b_spline = (t, knots, ind, order) => {\n const build_spline = (t, knots, ind, order) => {\n const b_spline = [];\n\n if(order == 1) {\n for(let i = 0; i < t.length; i++) {\n b_spline[i] = (knots[ind] <= t[i] && t[i] < knots[ind + 1]) ? 1 : 0;\n }\n }\n else {\n const b1 = build_spline(t, knots, ind, order - 1);\n const b2 = build_spline(t, knots, ind + 1, order - 1);\n\n for(let i = 0; i < t.length; i++) {\n let w1 = 0;\n let w2 = 0;\n if(knots[ind] != knots[ind + order - 1]) {\n w1 = (t[i] - knots[ind]) / (knots[ind + order - 1] - knots[ind]);\n }\n if(knots[ind + 1] != knots[ind + order]) {\n w2 = (knots[ind + order] - t[i]) / (knots[ind + order] - knots[ind + 1]);\n }\n\n b_spline[i] = w1 * b1[i] + w2 * b2[i];\n }\n }\n return b_spline;\n }\n return build_spline(t, knots, ind, order);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":261,"value":"indexOfMax = (arr) => arr.reduce((iMax, x, i, arr) => x > arr[iMax] ? i : iMax, 0);","pinned":false,"hidden":false,"mode":"ojs"},{"id":17,"value":"d3 = import(\"d3@v7\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":50,"value":"mathjs = require('mathjs@7.6.0')","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W37","stdlib":"1","titling":"auto"}