Code I wrote that I didn't happen to need for now

linear_regression.js

/*
 *  EOLIAN Web Music Player
 *  Copyright (C) 2018  Andreas Stöckel
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Affero General Public License as
 *  published by the Free Software Foundation, either version 3 of the
 *  License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Affero General Public License for more details.
 *
 *  You should have received a copy of the GNU Affero General Public License
 *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 */

/**
 * @file script/utils/linear_regression.js
 *
 * Implements simple linear regression. Used in conjunction with network time
 * synchronisation.
 *
 * @author Andreas Stöckel
 */

this.eolian = this.eolian || {};
this.eolian.utils = this.eolian.utils || {};
this.eolian.utils.linear_regression = (function (global) {
	'use strict';

	/**
	 * For the given x, y pairs computes the parameters of a line mx + b that
	 * minimizes quadratic error. This problem can be phrased as follows:
	 *
	 * y⃑ = A w⃑
	 *
	 * where y⃑ is the vector with all y-coordinates, A is a n x 2 matrix
	 * containing ones in the first column (for the bias term) and the vector
	 * of all x-coordinates in the second column, and w⃑ is the 2 x 1 vector
	 * (b, m).
	 *
	 * Now, the regularised least squares solution for this problem is
	 *
	 * w⃑ = (A^T A + λI)^-1 A^T y⃑
	 *
	 * The matrix A^T A is a 2 x 2 matrix, so the inverse is trivial to compute.
	 *
	 * @param xs: array of x-values.
	 * @param ys: array of y-values.
	 * @return b, m, rmse
	 */
	function linear_regression(xs, ys) {
		// Compute A^T A, track the maximum and minimum value in xs
		if ((!xs.length) || (!ys.length) || (xs.length !== ys.length)) {
			throw "xs and ys must be arrays and have the same non-zero length";
		}
		let n = xs.length;
		let a11 = n + 0.0, a12 = 0.0, a22 = 0.0; // a21 = a12
		for (let i = 0; i < n; i++) {
			a12 += xs[i] + 0.0;
			a22 += xs[i] * xs[i] + 0.0;
		}

		// Invert A
		let det = a11 * a22 - a12 * a12;
		let ai11 = a22 / det, ai22 = a11 / det, ai12 = -a12 / det;

		// Compute A^T y⃑
		let aty1 = 0.0, aty2 = 0.0;
		for (let i = 0; i < n; i++) {
			aty1 += ys[i];
			aty2 += ys[i] * xs[i];
		}

		// Compute w⃑ = (A^T A + λI)^-1 A^T y⃑
		const b = ai11 * aty1 + ai12 * aty2;
		const m = ai12 * aty1 + ai22 * aty2;

		// Compute the RMSE
		let rmse = 0.0;
		for (let i = 0; i < n; i++) {
			const e = (m * xs[i] + b - ys[i]);
			rmse += e * e;
		}
		rmse = Math.sqrt(rmse / n);

		return [b, m, rmse];
	}

	return linear_regression;
})(this);