Skip to content

Instantly share code, notes, and snippets.

@jonahwilliams

jonahwilliams/LinearRegression.js

Last active August 19, 2021 13:35
Show Gist options
  • Save jonahwilliams/62be9996afe5c2531625 to your computer and use it in GitHub Desktop.
Save jonahwilliams/62be9996afe5c2531625 to your computer and use it in GitHub Desktop.
Download ZIP
Smoothing Splines
Raw
README.md

Smoothing Splines avoid the knot selection choices of other spline methods by using the maximum number of knots, one per data point. Using a ridge regression penalty controlls overfitting and simplifies hyper-parameter selection to lambda alone.

Raw
index.html
<!DOCTYPE html>
<meta charset="utf-8">
<head>
<style>
.axis {
font: 10px sans-serif;
}
path {
stroke: steelblue;
stroke-width: 2;
fill: none;
}
.errors {
fill: red;
fill-opacity: 0.4;
}
.axis path,
.axis line {
fill: none;
stroke: #000;
shape-rendering: crispEdges;
}
</style>
</head>
<body>
<form>
<label id="label">lambda: 0.00001</label>
<input type="range" id="points" min="0" max="50" value="0" step="1">
</form>
<p id="dof">Estimated Degrees of Freedom: </p>
<script src="LinearRegression.js"></script>
<script src="Splines.js"></script>
<script src="http://d3js.org/d3.v3.min.js"></script>
<script>
var data = [];
for (var i = 0; i < 20; i++){
var a = Math.round(600 * Math.random())/100;
data.push({'x': a,'y': 0.9 * Math.sin(a) + (0.25 * Math.random()) - 0.125});
}
var lambda = 0.00001 * Math.pow(2, +document.getElementById("points").value);
var basis = smoothingSplines(data, data);
var output = RidgeRegression(basis, lambda),
params = output[0],
dof = output[1];
dof = Math.round(dof * 100)/100;
d3.select("#dof")
.text("Estimated Degrees of Freedom: " + dof);
// Generate smooth line
var Xhat = [],
i = 0,
j = 0;
for (i = 1; i < 600; i++){
var val = smoothingSplines([{'x': i/100, 'y':0}], data);
var yhat = 0.0;
for (j = 0; j < params.length; j++){
yhat += params[j][0] * val.x[0][j];
}
Xhat.push({'x': i/100, 'y': yhat});
}
var margin = {top: 80, right: 180, bottom: 80, left: 180},
width = 960 - margin.left - margin.right,
height = 500 - margin.top - margin.bottom;
var svg = d3.select("body").append("svg")
.attr("width", width + margin.left + margin.right)
.attr("height", height + margin.top + margin.bottom)
.append("g")
.attr("transform", "translate(" + margin.left + "," + margin.top + ")");
var y = d3.scale.linear()
.domain([-1, 1])
.range([height, 0]);
var x = d3.scale.linear()
.domain([0, 6])
.range([0, width])
var xAxis = d3.svg.axis()
.scale(x)
.orient("bottom");
var yAxis = d3.svg.axis()
.scale(y)
.orient("left");
var line = d3.svg.line()
.x(function(d) { return x(d.x); })
.y(function(d) { return y(d.y); });
svg.append("g")
.attr("class", "x axis")
.attr("transform", "translate(0," + height + ")")
.call(xAxis)
svg.append("g")
.attr("class", "y axis")
.call(yAxis);
svg.selectAll(".dot")
.data(data)
.enter().append("circle")
.attr("class", "dot")
.attr("r", 3.5)
.attr("cx", function(d){
return x(d.x);
})
.attr("cy", function(d){
return y(d.y)
})
.style("fill", "black");
svg.append("path")
.datum(Xhat)
.attr("class", "line")
.attr("d", line);
/* svg.selectAll("errors")
.data(errors)
.enter().append("rect")
.attr("class", "errors")
.attr("x", function(d){
return x(d.x);
})
.attr("y", function(d){
if (d.sy <= d.y){
return y(d.y);
}
else {
return y(d.y) - Math.abs(y(d.y) - y(d.sy));
}
})
.attr("width", function(d){
return Math.abs(x(d.y) - x(d.sy));
})
.attr("height", function(d){
return Math.abs(y(d.y) - y(d.sy));
})
.style("fill", "red")
.style("fill-opacity", 0.5);
*/
d3.select("#points")
.on("change", function(d){
var lambda = 0.00001 * Math.pow(2, +document.getElementById("points").value);
d3.select("label")
.text("lambda: " + lambda);
var output = RidgeRegression(basis, +lambda),
params = output[0],
dof = output[1];
dof = Math.round(dof * 100)/100;
d3.select("#dof")
.text("Estimated Degrees of Freedom: " + dof);
var Xhat = [],
i = 0,
j = 0;
for (i = 1; i < 600; i++){
var val = smoothingSplines([{'x': i/100, 'y':0}], data);
var yhat = 0.0;
for (j = 0; j < params.length; j++){
yhat += params[j][0] * val.x[0][j];
}
Xhat.push({'x': i/100, 'y': yhat});
}
svg.selectAll(".line")
.datum(Xhat)
.transition()
.attr("class", "line")
.attr("d", line);
/* svg.selectAll(".errors")
.data(errors)
.transition()
.attr("class","errors")
.attr("x", function(d){
return x(d.x);
})
.attr("y", function(d){
if (d.sy <= d.y){
return y(d.y);
}
else {
return y(d.y) - Math.abs(y(d.y) - y(d.sy));
}
})
.attr("width", function(d){
return Math.abs(x(d.y) - x(d.sy));
})
.attr("height", function(d){
return Math.abs(y(d.y) - y(d.sy));
})
.style("fill", "red")
.style("fill-opacity", 0.5);
*/
});
</script>
</body>
Raw
LinearRegression.js
function LinearRegression(data){
var X = [],
y = [];
// Get everything in array format
for (var i = 0; i < data.length; i ++){
X[i] = [];
y[i] = data[i].y;
for(var j = 0; j < data[0]['x'].length; j++){
X[i][j] = data[i]['x'][j];
}
}
var t = matrixTranspose(X);
var hatMatrix = matrixMultiply(matrixMultiply(matrixInvert(matrixMultiply(t, X)), t), y);
return hatMatrix;
}
function RidgeRegression(data, lambda){
var X = data.x,
y = data.y;
var size = X[0].length,
t = matrixTranspose(X),
reg = matrixRegularizer(size, lambda)
var inner = matrixAdd(matrixMultiply(t, X), reg);
var hatMatrix = matrixMultiply(matrixMultiply(matrixInvert(inner), t), y);
var dof = matrixTrace(matrixMultiply(X, matrixMultiply(matrixInvert(inner), t)));
return [hatMatrix, dof];
}
function matrixTranspose(a){
//http://stackoverflow.com/questions/4492678/to-swap-rows-with-columns-of-matrix-in-javascript-or-jquery
return Object.keys(a[0]).map(
function (c) { return a.map(function (r) { return r[c]; }); }
);
}
function matrixMultiply(m1, m2) {
// http://tech.pro/tutorial/1527/matrix-multiplication-in-functional-javascript
//
var result = [];
for (var i = 0; i < m1.length; i++) {
result[i] = [];
for (var j = 0; j < m2[0].length; j++) {
var sum = 0;
for (var k = 0; k < m1[0].length; k++) {
sum += m1[i][k] * m2[k][j];
}
result[i][j] = sum;
}
}
return result;
}
function matrixInvert(M){
// source - http://blog.acipo.com/matrix-inversion-in-javascript/
//
if(M.length !== M[0].length){return;}
var i=0, ii=0, j=0, dim=M.length, e=0, t=0;
var I = [], C = [];
for(i=0; i<dim; i+=1){
I[I.length]=[];
C[C.length]=[];
for(j=0; j<dim; j+=1){
if(i==j){ I[i][j] = 1; }
else{ I[i][j] = 0; }
C[i][j] = M[i][j];
}
}
for(i=0; i<dim; i+=1){
e = C[i][i];
if(e==0){
for(ii=i+1; ii<dim; ii+=1){
if(C[ii][i] != 0){
for(j=0; j<dim; j++){
e = C[i][j];
C[i][j] = C[ii][j];
C[ii][j] = e;
e = I[i][j];
I[i][j] = I[ii][j];
I[ii][j] = e;
}
break;
}
}
e = C[i][i];
if(e==0){return}
}
for(j=0; j<dim; j++){
C[i][j] = C[i][j]/e;
I[i][j] = I[i][j]/e;
}
for(ii=0; ii<dim; ii++){
if(ii==i){continue;}
e = C[ii][i];
for(j=0; j<dim; j++){
C[ii][j] -= e*C[i][j];
I[ii][j] -= e*I[i][j];
}
}
}
return I;
}
function matrixRegularizer(size, lambda){
//Returns the identity matrix for a given size times lambda:
// First generate a zero filled matrix
var regularizer = [];
var i = 0, j = 0;
for(i = 0; i < size; i++){
regularizer[i] = [];
for (j = 0; j < size; j++){
regularizer[i][j] = 0;
}
}
//Then populate the diagonal with lambda
for (i = 0; i < size; i++){
regularizer[i][i] = lambda;
}
return regularizer;
}
function matrixAdd(m1, m2){
//Adds two matrixes of the same size. DOES NOT CHECK FOR THIS
//uses matrix 1 as the new matrix
for (var i = 0; i < m1.length; i++){
for (var j = 0; j < m1[0].length; j++){
m1[i][j] = m1[i][j] + m2[i][j];
}
};
return m1;
}
function matrixTrace(m){
// Calculates the trace of the matrix for estimating degrees of freedom
// Assumes matrix is square
var i = 0,
j = 0,
n = m.length,
trace = 0.0;
for (i = 0; i < n; i++){
trace += m[i][i];
}
return trace;
}
Raw
Splines.js
function smoothingSplines(data, spline_locations){
var i = 0,
j = 0;
var X = [],
y = [],
n = data.length,
p = spline_locations.length;
//Allocate Inital Matrix
for (i = 0; i < n; i++){
X[i] = [];
for (j = 0; j < p + 4; j++){
X[i][j] = 0;
}
}
for (i = 0; i < n; i++){
X[i][0] = 1;
X[i][1] = data[i].x;
X[i][2] = Math.pow(data[i].x, 2);
X[i][3] = Math.pow(data[i].x, 3);
y[i] = [data[i].y];
//Now allocate knot locations
for (j = 0; j < p; j++){
X[i][j + 3] = Math.max(
Math.pow(X[i][1] - spline_locations[j].x, 3), 0
);
}
}
return {'x': X, 'y': y};
}
@syntheticzero
Copy link

Is this code available for free use? I'm interested in adapting this for my project. BSD License? MIT? I guess I'm assuming it is, since it's posted here, but if you have any specific comments about licensing please let me know (comment here). Thanks...

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment