31:T949,function fitLeastSquares(x, y, basisFns) { if (x.length !== y.length) { throw new Error("Incompatible input vectors, x.length != y.length"); } if (x.length < basisFns.length) { throw new Error("Insufficient data, x.length < basisFns.length"); } const numBasis = basisFns.length; const N = x.length; // Design matrix A (row-major flat array of size N x M) const A = Array(N * numBasis).fill(0); for (let j = 0; j < numBasis; j++) { for (let i = 0; i < N; i++) { A[i * numBasis + j] = basisFns[j](x[i]); } } function column(A, rows, cols, j) { const col = Array(rows).fill(0); for (let i = 0; i < rows; i++) col[i] = A[i * cols + j]; return col; } function dot(a, b) { let sum = 0; for (let i = 0; i < a.length; i++) sum += a[i] * b[i]; return sum; } // QR Decomposition (Gram-Schmidt) function gramSchmidtQRDecomp(A, N, M) { const Q = Array(N * M).fill(0); const R = Array(M * M).fill(0); for (let j = 0; j < M; j++) { const v = column(A, N, M, j); for (let i = 0; i < j; i++) { const qi = column(Q, N, M, i); const rij = dot(qi, v); R[i * M + j] = rij; for (let k = 0; k < N; k++) { v[k] -= rij * qi[k]; } } const rjj = Math.sqrt(dot(v, v)); R[j * (M + 1)] = rjj; for (let k = 0; k < N; k++) { Q[k * M + j] = v[k] / rjj; } } return { Q, R }; } // Back substitution function backSubstitute(R, Qt_y, M) { const beta = Array(M).fill(0); for (let i = M - 1; i >= 0; i--) { let sum = Qt_y[i]; for (let j = i + 1; j < M; j++) { sum -= R[i * M + j] * beta[j]; } beta[i] = sum / R[i * (M + 1)]; } return beta; } const { Q, R } = gramSchmidtQRDecomp(A, N, numBasis); // Compute Q^T y const Qt_y = Array(numBasis).fill(0); for (let i = 0; i < N; i++) { for (let j = 0; j < numBasis; j++) { Qt_y[j] += Q[i * numBasis + j] * y[i]; } } const coefficients = backSubstitute(R, Qt_y, numBasis); let sumSquares = 0.0; for (let i = 0; i < N; i++) { let Ax_i = 0; for (let j = 0; j < numBasis; j++) { Ax_i += A[i * numBasis + j] * coefficients[j]; } const residual = y[i] - Ax_i; sumSquares += residual ** 2; } return { coefficients, error: Math.sqrt(sumSquares / N) }; }32:T609,signature(fitLeastSquares, { description: md` Use a [QR decomposition](https://en.wikipedia.org/wiki/QR_decomposition) to perform a least squares fit of arbitrary basis functions to a set of ${tex`(x_i, y_i)`} input data. **Import**: \`\`\`js import { fitLeastSquares } from '${await getPinnedSlug( "@rreusser/least-squares-curve-fit" )}' \`\`\` **Parameters**: - \`x: number[]\`: ${tex`x`}-coordinate of input data - \`y: number[]\`: ${tex`y`}-coordinate of input data - \`basisFns: Array<(number) => number>\`: list of basis functions **Returns**: - \`coefficients: number[]\`: list of coefficients for the respective basis functions. - \`error: number\`: root mean square error (RMSE) **Notes**: - The length of arrays \`x\` and \`y\` must be equal. - Length of input data must be greater than or equal to the number of basis functions.`, example: [ ` // Fit a quadratic polynomial to y = 1 + x² fitLeastSquares( [0, 1, 2, 3, 4], [1, 2, 5, 10, 17], [(x) => 1, (x) => x, (x) => x**2] ) // => { coefficients: [ 1, 0, 1 ], error: 0 }` ], tests: { "basic curve fit": (assert) => { const fit = fitLeastSquares( [0, 1, 2, 3, 4], [1, 2, 5, 10, 17], [(x) => 1, (x) => x, (x) => x ** 2] ); const expected = [1, 0, 1]; assert(Math.abs(fit.coefficients[0] - expected[0]) < 1e-8); assert(Math.abs(fit.coefficients[1] - expected[1]) < 1e-8); assert(Math.abs(fit.coefficients[2] - expected[2]) < 1e-8); assert(Math.abs(fit.error - 0 < 1e-8)); } } })33:[] 2f:{"title":"Least Squares Curve Fit with the QR Decomposition","cells":[{"id":0,"value":"# Least Squares Curve Fit with the QR Decomposition","pinned":false,"hidden":false,"mode":"md"},{"id":4,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":219,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":281,"value":"## Example: Fit a fifth-order polynomial to ${tex`\\sin(x)`}","pinned":false,"hidden":false,"mode":"md"},{"id":302,"value":"N = 51","pinned":true,"hidden":false,"mode":"ojs"},{"id":45,"value":"x = Array.from({ length: N }, (_, i) => (Math.PI * 2 * i) / (N - 1))","pinned":true,"hidden":false,"mode":"ojs"},{"id":87,"value":"y = x.map(x => Math.sin(x))","pinned":true,"hidden":false,"mode":"ojs"},{"id":78,"value":"basisFns = [\n (x) => x ** 0,\n (x) => x ** 1,\n (x) => x ** 2,\n (x) => x ** 3,\n (x) => x ** 4,\n (x) => Math.exp(x) // Doesn't have to be a polynomial!\n]","pinned":true,"hidden":false,"mode":"ojs"},{"id":60,"value":"fit = fitLeastSquares(x, y, basisFns)","pinned":true,"hidden":false,"mode":"ojs"},{"id":82,"value":"{\n const N = 101;\n const yFit = x.map((x) =>\n basisFns\n .slice(0, fit.coefficients.length)\n .reduce((sum, fn, i) => sum + fit.coefficients[i] * fn(x), 0)\n );\n const plotData = x.map((x, i) => ({\n x: x,\n y: y[i],\n yFit: yFit[i]\n }));\n\n return Plot.plot({\n grid: true,\n marks: [\n Plot.dot(plotData, {\n x: \"x\",\n y: \"y\",\n r: 2.5,\n fill: \"black\",\n stroke: \"none\"\n }),\n Plot.line(plotData, {\n x: \"x\",\n y: \"yFit\",\n stroke: \"red\",\n strokeWidth: 2\n })\n ],\n marginLeft: 40,\n marginBottom: 40\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":215,"value":"import { signature, getPinnedSlug } from \"@mootari/signature@553\"","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W33","stdlib":"1","titling":"auto"}