Skip to content

Instantly share code, notes, and snippets.

@armollica
Last active April 15, 2016 02:16
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save armollica/56217d01ddf1370773da to your computer and use it in GitHub Desktop.
Save armollica/56217d01ddf1370773da to your computer and use it in GitHub Desktop.
Multivariate linear regression
height: 300

Click on the chart to change the model prediction data point (orange line) and see how the independent variables impact one another's marginal effect on the dependent variable.

Shown are the model predictions (solid line) and 95% confidence regression interval (dashed line) for the following model:

mpg = B0 + B1 wt + B2 hp + B3 (wt x hp) + B4 wt2 + B5 hp2

where

mph = Miles/Gallon

wt = Weight (1000 lb)

hp = Gross horsepower

Because this model is linear in the parameters, Bj, it can be estimated using least-squares. The independent variables, weight and horsepower, enter the model nonlinearly and so the partial derivative of weight and horsepower will vary as weight and horsepower change. In other words, the "marginal effect" on fuel efficiency from an increase in a car's weight is different depending upon whether the car was intially light or heavy.

The slopes on the above graphs show this marginal effect. The downward sloping functions confirm what would be expected; cars that weigh more and have more horsepower are predicted to be less fuel efficient (have lower miles/gallon).

The (slight) nonlinearity in the regression equation shows that an increase in a car's horsepower has a varying impact on fuel efficiency depending on whether the car initially had a lot of horsepower or just a little. The slope is steeper at lower horsepower levels and flatter at higher horsepower levels. This implies that an increasing a car's horsepower by 20 is predicted to decrease fuel efficiency more for a low-horsepower car than it is for a high-horsepower car.

Because of the cross-product term, (wt x hp), the marginal effect of a car's weight on miles/gallon depends upon the horsepower of the car and vice versa. Clicking on the chart changes the weight and horsepower. Changing the weight changes the marginal effect for horsepower and changing the horsepower changes the marginal effect for weight. Increasing a car's weight makes the regression equation for horsepower flatter and lower. This implies that for heavier cars, horsepower plays a smaller role in determining fuel efficiency. Note that the orange indicator line cross the model prediction line at the same place.

Here are the equations for the marginal effects (i.e., partial derivatives of mpg):

mpg'wt = B1 + B3 hp + 2 B4 wt

mpg'hp = B2 + B3 wt + 2 B5 hp

Notes: The model and data used here are for demonstration purposes and the model predictions are very weak approximations to a complex system. To see the model results open up the console log in the browser (F12 on Chrome). The data source is mtcars from R's dataset package. Formulas for regression interval and White's scaled variance-covariance matrix were found here: http://www.ssc.wisc.edu/~bhansen/econometrics/ Using science.js for linear algebra functions.

[{"mpg":21,"cyl":6,"disp":160,"hp":110,"drat":3.9,"wt":2.62,"qsec":16.46,"vs":0,"am":1,"gear":4,"carb":4,"_row":"Mazda RX4"},{"mpg":21,"cyl":6,"disp":160,"hp":110,"drat":3.9,"wt":2.875,"qsec":17.02,"vs":0,"am":1,"gear":4,"carb":4,"_row":"Mazda RX4 Wag"},{"mpg":22.8,"cyl":4,"disp":108,"hp":93,"drat":3.85,"wt":2.32,"qsec":18.61,"vs":1,"am":1,"gear":4,"carb":1,"_row":"Datsun 710"},{"mpg":21.4,"cyl":6,"disp":258,"hp":110,"drat":3.08,"wt":3.215,"qsec":19.44,"vs":1,"am":0,"gear":3,"carb":1,"_row":"Hornet 4 Drive"},{"mpg":18.7,"cyl":8,"disp":360,"hp":175,"drat":3.15,"wt":3.44,"qsec":17.02,"vs":0,"am":0,"gear":3,"carb":2,"_row":"Hornet Sportabout"},{"mpg":18.1,"cyl":6,"disp":225,"hp":105,"drat":2.76,"wt":3.46,"qsec":20.22,"vs":1,"am":0,"gear":3,"carb":1,"_row":"Valiant"},{"mpg":14.3,"cyl":8,"disp":360,"hp":245,"drat":3.21,"wt":3.57,"qsec":15.84,"vs":0,"am":0,"gear":3,"carb":4,"_row":"Duster 360"},{"mpg":24.4,"cyl":4,"disp":146.7,"hp":62,"drat":3.69,"wt":3.19,"qsec":20,"vs":1,"am":0,"gear":4,"carb":2,"_row":"Merc 240D"},{"mpg":22.8,"cyl":4,"disp":140.8,"hp":95,"drat":3.92,"wt":3.15,"qsec":22.9,"vs":1,"am":0,"gear":4,"carb":2,"_row":"Merc 230"},{"mpg":19.2,"cyl":6,"disp":167.6,"hp":123,"drat":3.92,"wt":3.44,"qsec":18.3,"vs":1,"am":0,"gear":4,"carb":4,"_row":"Merc 280"},{"mpg":17.8,"cyl":6,"disp":167.6,"hp":123,"drat":3.92,"wt":3.44,"qsec":18.9,"vs":1,"am":0,"gear":4,"carb":4,"_row":"Merc 280C"},{"mpg":16.4,"cyl":8,"disp":275.8,"hp":180,"drat":3.07,"wt":4.07,"qsec":17.4,"vs":0,"am":0,"gear":3,"carb":3,"_row":"Merc 450SE"},{"mpg":17.3,"cyl":8,"disp":275.8,"hp":180,"drat":3.07,"wt":3.73,"qsec":17.6,"vs":0,"am":0,"gear":3,"carb":3,"_row":"Merc 450SL"},{"mpg":15.2,"cyl":8,"disp":275.8,"hp":180,"drat":3.07,"wt":3.78,"qsec":18,"vs":0,"am":0,"gear":3,"carb":3,"_row":"Merc 450SLC"},{"mpg":10.4,"cyl":8,"disp":472,"hp":205,"drat":2.93,"wt":5.25,"qsec":17.98,"vs":0,"am":0,"gear":3,"carb":4,"_row":"Cadillac Fleetwood"},{"mpg":10.4,"cyl":8,"disp":460,"hp":215,"drat":3,"wt":5.424,"qsec":17.82,"vs":0,"am":0,"gear":3,"carb":4,"_row":"Lincoln Continental"},{"mpg":14.7,"cyl":8,"disp":440,"hp":230,"drat":3.23,"wt":5.345,"qsec":17.42,"vs":0,"am":0,"gear":3,"carb":4,"_row":"Chrysler Imperial"},{"mpg":32.4,"cyl":4,"disp":78.7,"hp":66,"drat":4.08,"wt":2.2,"qsec":19.47,"vs":1,"am":1,"gear":4,"carb":1,"_row":"Fiat 128"},{"mpg":30.4,"cyl":4,"disp":75.7,"hp":52,"drat":4.93,"wt":1.615,"qsec":18.52,"vs":1,"am":1,"gear":4,"carb":2,"_row":"Honda Civic"},{"mpg":33.9,"cyl":4,"disp":71.1,"hp":65,"drat":4.22,"wt":1.835,"qsec":19.9,"vs":1,"am":1,"gear":4,"carb":1,"_row":"Toyota Corolla"},{"mpg":21.5,"cyl":4,"disp":120.1,"hp":97,"drat":3.7,"wt":2.465,"qsec":20.01,"vs":1,"am":0,"gear":3,"carb":1,"_row":"Toyota Corona"},{"mpg":15.5,"cyl":8,"disp":318,"hp":150,"drat":2.76,"wt":3.52,"qsec":16.87,"vs":0,"am":0,"gear":3,"carb":2,"_row":"Dodge Challenger"},{"mpg":15.2,"cyl":8,"disp":304,"hp":150,"drat":3.15,"wt":3.435,"qsec":17.3,"vs":0,"am":0,"gear":3,"carb":2,"_row":"AMC Javelin"},{"mpg":13.3,"cyl":8,"disp":350,"hp":245,"drat":3.73,"wt":3.84,"qsec":15.41,"vs":0,"am":0,"gear":3,"carb":4,"_row":"Camaro Z28"},{"mpg":19.2,"cyl":8,"disp":400,"hp":175,"drat":3.08,"wt":3.845,"qsec":17.05,"vs":0,"am":0,"gear":3,"carb":2,"_row":"Pontiac Firebird"},{"mpg":27.3,"cyl":4,"disp":79,"hp":66,"drat":4.08,"wt":1.935,"qsec":18.9,"vs":1,"am":1,"gear":4,"carb":1,"_row":"Fiat X1-9"},{"mpg":26,"cyl":4,"disp":120.3,"hp":91,"drat":4.43,"wt":2.14,"qsec":16.7,"vs":0,"am":1,"gear":5,"carb":2,"_row":"Porsche 914-2"},{"mpg":30.4,"cyl":4,"disp":95.1,"hp":113,"drat":3.77,"wt":1.513,"qsec":16.9,"vs":1,"am":1,"gear":5,"carb":2,"_row":"Lotus Europa"},{"mpg":15.8,"cyl":8,"disp":351,"hp":264,"drat":4.22,"wt":3.17,"qsec":14.5,"vs":0,"am":1,"gear":5,"carb":4,"_row":"Ford Pantera L"},{"mpg":19.7,"cyl":6,"disp":145,"hp":175,"drat":3.62,"wt":2.77,"qsec":15.5,"vs":0,"am":1,"gear":5,"carb":6,"_row":"Ferrari Dino"},{"mpg":15,"cyl":8,"disp":301,"hp":335,"drat":3.54,"wt":3.57,"qsec":14.6,"vs":0,"am":1,"gear":5,"carb":8,"_row":"Maserati Bora"},{"mpg":21.4,"cyl":4,"disp":121,"hp":109,"drat":4.11,"wt":2.78,"qsec":18.6,"vs":1,"am":1,"gear":4,"carb":2,"_row":"Volvo 142E"}]
<html>
<head>
<style>
body {
font: 12px sans-serif;
}
.axis path,
.axis line {
fill: none;
stroke: black;
stroke-width: 1px;
}
.line {
fill: none;
stroke: black;
}
.line.upper,
.line.lower {
stroke-dasharray: 5, 5;
}
.marker {
stroke: tomato;
stroke-width: 2px;
}
</style>
</head>
<body>
<script src="//d3js.org/d3.v3.min.js" charset="utf-8"></script>
<script src="science.v1.min.js"></script>
<script src="lm.js"></script>
<script src="line-chart.js"></script>
<script>
var margin = { top: 10, left: 30, bottom: 50, right: 10 },
width = 480 - margin.left - margin.right,
height = 300 - margin.top - margin.bottom;
var lineChart = d3.lineChart() // see line-chart.js
.width(width)
.height(height)
.yExtent([0, 32]);
[ { key: "wt", name: "Weight (1000 lbs)"},
{ key: "hp", name: "Gross horsepower" }]
.forEach(function(d) {
var container = d3.select("body").append("div")
.attr("class", "container " + d.key)
.style("position", "relative")
.style("float", "left");
var svg = container.append("svg")
.attr("width", width + margin.left + margin.right)
.attr("height", height + margin.top + margin.bottom)
.append("g")
.attr("class", "canvas")
.attr("transform", "translate(" + margin.left + "," + margin.top + ")");
svg.append("line")
.attr("class", "marker");
svg.append("g").attr("class", "x axis")
.attr("transform", "translate(0," + height + ")")
.append("text")
.attr("class", "label")
.attr("x", width)
.attr("y", -6)
.style("text-anchor", "end")
.text(d.name);
svg.append("g").attr("class", "y axis")
.append("text")
.attr("class", "label")
.attr("y", 6)
.attr("dx", ".71em")
.style("text-anchor", "start")
.text("Miles/Gallon")
svg.append("rect")
.attr("class", "event-rect")
.attr("width", width)
.attr("height", height)
.style("fill-opacity", 0);
});
d3.json("cars.json", function(error, cars) {
if (error) throw error;
var model = lm() // see lm.js
.y(function(d) { return d.mpg; })
.X(function(d) {
return [1, d.wt, d.hp, d.wt*d.hp, Math.pow(d.wt, 2), Math.pow(d.hp, 2)];
})
.data(cars);
var fit = model();
console.log(fit);
var wt = d3.mean(cars, pluck("wt"));
var wtExtent = d3.extent(cars, pluck("wt"));
var hp = d3.mean(cars, pluck("hp"));
var hpExtent = d3.extent(cars, pluck("hp"));
d3.select(".wt .event-rect")
.on("click", function() {
var x = d3.mouse(this)[0];
var wtInterval = renderWtInterval(hp);
wt = wtInterval.xScale.invert(x);
var hpInterval = renderHpInterval(wt);
d3.select(".wt .marker")
.attr("x1", x)
.attr("x2", x);
});
d3.select(".hp .event-rect")
.on("click", function() {
var x = d3.mouse(this)[0];
var hpInterval = renderHpInterval(wt);
hp = hpInterval.xScale.invert(x);
var wtInterval = renderWtInterval(hp);
d3.select(".hp .marker")
.attr("x1", x)
.attr("x2", x);
});
var wtInterval = renderWtInterval(hp);
var wtScale = wtInterval.xScale;
d3.select(".wt .marker")
.attr("x1", wtScale(wt))
.attr("x2", wtScale(wt))
.attr("y1", height)
.attr("y2", 0);
var hpInterval = renderHpInterval(wt);
var hpScale = hpInterval.xScale;
d3.select(".hp .marker")
.attr("x1", hpScale(hp))
.attr("x2", hpScale(hp))
.attr("y1", height)
.attr("y2", 0);
function renderWtInterval(hp) {
var wt_data = sequence(cars, pluck("wt"), 40)
.map(function(wt) {
var interval = fit.getRegressionInterval({ wt: wt, hp: hp });
return {
wt: wt,
y_pred: interval.y_pred,
y_upper: interval.y_upper,
y_lower: interval.y_lower
};
});
wtChart = lineChart
.x(function(d) { return d.wt; });
d3.select(".wt svg").select(".canvas")
.call(wtChart.y(pluck("y_pred")), "y line", wt_data)
.call(wtChart.y(pluck("y_lower")), "lower line", wt_data)
.call(wtChart.y(pluck("y_upper")), "upper line", wt_data);
d3.select(".wt svg").select(".x.axis").call(wtChart.xAxis());
d3.select(".wt svg").select(".y.axis").call(wtChart.yAxis());
return {
data: wt_data,
xScale: wtChart.xScale(),
yScale: wtChart.yScale()
};
}
function renderHpInterval(wt) {
var wt_mean = d3.mean(cars, pluck("wt"));
var hp_data = sequence(cars, pluck("hp"), 40)
.map(function(hp) {
var interval = fit.getRegressionInterval({ wt: wt, hp: hp });
return {
hp: hp,
y_pred: interval.y_pred,
y_upper: interval.y_upper,
y_lower: interval.y_lower
};
});
hpChart = lineChart
.x(function(d) { return d.hp; });
d3.select(".hp svg").select(".canvas")
.call(hpChart.y(pluck("y_pred")), "y line", hp_data)
.call(hpChart.y(pluck("y_lower")), "lower line", hp_data)
.call(hpChart.y(pluck("y_upper")), "upper line", hp_data);
d3.select(".hp svg").select(".x.axis").call(hpChart.xAxis());
d3.select(".hp svg").select(".y.axis").call(hpChart.yAxis());
return {
data: hp_data,
xScale: hpChart.xScale(),
yScale: hpChart.yScale()
};
}
});
// Returns accessor function for a given key
// makes accessing "columns" from associative arrays easier
function pluck(key) { return function(d) { return d[key]; }; }
// Returns a sequence of numbers (of specified length) between the min
// and max of an array
function sequence(data, accessor, length_out) {
var extent = d3.extent(data, accessor);
var step = (extent[1] - extent[0])/length_out;
return d3.range(extent[0], extent[1], step);
}
</script>
</body>
</html>
d3.lineChart = function() {
var width = 400,
height = 300,
xExtent, yExtent,
x, y,
xScale = d3.scale.linear(),
yScale = d3.scale.linear(),
path = d3.svg.line();
function chart(selection, lineClass, data) {
xScale
.domain(xExtent || d3.extent(data, x))
.range([0, width]);
yScale
.domain(yExtent || d3.extent(data, y))
.range([height, 0]);
path
.x(function(d) { return xScale(x(d)); })
.y(function(d) { return yScale(y(d)); });
var lineSelectorString = "." + (lineClass.replace(" ", "."));
var lines = selection.selectAll(lineSelectorString).data([data]);
lines.enter().append("path")
.attr("class", lineClass);
lines
.transition().duration(333)
.attr("d", path);
lines.exit().remove();
}
chart.width = function(_) {
if (!arguments.length) return width;
width = _;
return chart;
};
chart.height = function(_) {
if (!arguments.length) return height;
height = _;
return chart;
};
chart.x = function(_) {
if (!arguments.length) return x;
x = _;
return chart;
};
chart.y = function(_) {
if (!arguments.length) return y;
y = _;
return chart;
};
chart.xExtent = function(_) {
if (!arguments.length) return xExtent;
xExtent = _;
return chart;
};
chart.yExtent = function(_) {
if (!arguments.length) return yExtent;
yExtent = _;
return chart;
};
// Getter functions
chart.xAxis = function() { return d3.svg.axis().scale(xScale).orient("bottom"); };
chart.yAxis = function() { return d3.svg.axis().scale(yScale).orient("left"); };
chart.xScale = function() { return xScale; };
chart.yScale = function() { return yScale; };
return chart;
}
// Multivariate linear regression
function lm() {
//make science.lin function easier to reference
var m = science.lin.multiply,
t = science.lin.transpose,
inv = science.lin.inverse;
var _y, _X, _data;
function model() {
var y = vectorToMatrix(_data.map(_y)); // dependent variable (n x 1 matrix)
var X = _data.map(_X); // independent variables (n x k matrix)
var n = X.length; // n = # of observations
var k = X[0].length; // k = # of model parameters
var invXX = inv(m(t(X), X)); // (X'X)^1 (k x k matrix)
var B = m(invXX, m(t(X), y)); // model parameters (k x 1 matrix)
var y_pred = m(X, B); // in-sample model predictions
var e = d3.range(n).map(function(i) { // prediction error
return [y[i][0] - y_pred[i][0]];
});
// Variance-covariance matrix using scaled White's formula
// (see chapter 4 from http://www.ssc.wisc.edu/~bhansen/econometrics/)
var u = X.map(function(x, i) {
return x.map(function(d) { return d * e[i][0]; });
});
var whites_scaler = n/(n-k);
var V = m(invXX, m(m(t(u), u), invXX))
.map(function(row) {
return row.map(function(d) { return whites_scaler * d; });
});
// parameter standard errors
var se = diag(V).map(function(d) { return [Math.sqrt(d[0])]; });
// t-statistics
var tstat = B.map(function(b, i) { return [b[0] / se[i][0]]; });
return {
B: B,
se: se,
V: V,
t: tstat,
getRegressionInterval: getRegressionInterval
};
// get asymptotic 95% confidence regression interval
// at a given data point
function getRegressionInterval(d) {
var x = [_X(d)]; // 1 x k
var y_pred = m(x, B)[0][0];
var interval = 1.96 * Math.sqrt(m(m(x, V), t(x))[0][0]);
return {
y_lower: y_pred - interval,
y_pred: y_pred,
y_upper: y_pred + interval
};
}
}
model.y = function(_) {
if (!arguments.length) return _y;
_y = _;
return model;
};
model.X = function(_) {
if (!arguments.length) return _X;
_X = _;
return model;
};
model.data = function(_) {
if (!arguments.length) return _data;
_data = _;
return model;
};
return model;
// convert to/from "array of numbers" from/to "1D matrix"
function vectorToMatrix(_) { return _.map(function(d) { return [d]; }); }
function matrixToVector(_) { return _.map(function(d) { return d[0]; }); }
// create matrix of size r x c with all 1's
function ones(r, c) {
return d3.range(r).map(function() {
return d3.range(c).map(function() { return 1; });
});
}
// get diagonal vector from square matrix
function diag(x) { return x.map(function(d, j) { return [d[j]]; }); }
}
(function(exports){(function(exports){science={version:"1.9.3"};science.ascending=function(a,b){return a-b};science.EULER=.5772156649015329;science.expm1=function(x){return x<1e-5&&x>-1e-5?x+.5*x*x:Math.exp(x)-1};science.functor=function(v){return typeof v==="function"?v:function(){return v}};science.hypot=function(x,y){x=Math.abs(x);y=Math.abs(y);var max,min;if(x>y){max=x;min=y}else{max=y;min=x}var r=min/max;return max*Math.sqrt(1+r*r)};science.quadratic=function(){var complex=false;function quadratic(a,b,c){var d=b*b-4*a*c;if(d>0){d=Math.sqrt(d)/(2*a);return complex?[{r:-b-d,i:0},{r:-b+d,i:0}]:[-b-d,-b+d]}else if(d===0){d=-b/(2*a);return complex?[{r:d,i:0}]:[d]}else{if(complex){d=Math.sqrt(-d)/(2*a);return[{r:-b,i:-d},{r:-b,i:d}]}return[]}}quadratic.complex=function(x){if(!arguments.length)return complex;complex=x;return quadratic};return quadratic};science.zeroes=function(n){var i=-1,a=[];if(arguments.length===1)while(++i<n)a[i]=0;else while(++i<n)a[i]=science.zeroes.apply(this,Array.prototype.slice.call(arguments,1));return a}})(this);(function(exports){science.lin={};science.lin.decompose=function(){function decompose(A){var n=A.length,V=[],d=[],e=[];for(var i=0;i<n;i++){V[i]=[];d[i]=[];e[i]=[]}var symmetric=true;for(var j=0;j<n;j++){for(var i=0;i<n;i++){if(A[i][j]!==A[j][i]){symmetric=false;break}}}if(symmetric){for(var i=0;i<n;i++)V[i]=A[i].slice();science_lin_decomposeTred2(d,e,V);science_lin_decomposeTql2(d,e,V)}else{var H=[];for(var i=0;i<n;i++)H[i]=A[i].slice();science_lin_decomposeOrthes(H,V);science_lin_decomposeHqr2(d,e,H,V)}var D=[];for(var i=0;i<n;i++){var row=D[i]=[];for(var j=0;j<n;j++)row[j]=i===j?d[i]:0;D[i][e[i]>0?i+1:i-1]=e[i]}return{D:D,V:V}}return decompose};function science_lin_decomposeTred2(d,e,V){var n=V.length;for(var j=0;j<n;j++)d[j]=V[n-1][j];for(var i=n-1;i>0;i--){var scale=0,h=0;for(var k=0;k<i;k++)scale+=Math.abs(d[k]);if(scale===0){e[i]=d[i-1];for(var j=0;j<i;j++){d[j]=V[i-1][j];V[i][j]=0;V[j][i]=0}}else{for(var k=0;k<i;k++){d[k]/=scale;h+=d[k]*d[k]}var f=d[i-1];var g=Math.sqrt(h);if(f>0)g=-g;e[i]=scale*g;h=h-f*g;d[i-1]=f-g;for(var j=0;j<i;j++)e[j]=0;for(var j=0;j<i;j++){f=d[j];V[j][i]=f;g=e[j]+V[j][j]*f;for(var k=j+1;k<=i-1;k++){g+=V[k][j]*d[k];e[k]+=V[k][j]*f}e[j]=g}f=0;for(var j=0;j<i;j++){e[j]/=h;f+=e[j]*d[j]}var hh=f/(h+h);for(var j=0;j<i;j++)e[j]-=hh*d[j];for(var j=0;j<i;j++){f=d[j];g=e[j];for(var k=j;k<=i-1;k++)V[k][j]-=f*e[k]+g*d[k];d[j]=V[i-1][j];V[i][j]=0}}d[i]=h}for(var i=0;i<n-1;i++){V[n-1][i]=V[i][i];V[i][i]=1;var h=d[i+1];if(h!=0){for(var k=0;k<=i;k++)d[k]=V[k][i+1]/h;for(var j=0;j<=i;j++){var g=0;for(var k=0;k<=i;k++)g+=V[k][i+1]*V[k][j];for(var k=0;k<=i;k++)V[k][j]-=g*d[k]}}for(var k=0;k<=i;k++)V[k][i+1]=0}for(var j=0;j<n;j++){d[j]=V[n-1][j];V[n-1][j]=0}V[n-1][n-1]=1;e[0]=0}function science_lin_decomposeTql2(d,e,V){var n=V.length;for(var i=1;i<n;i++)e[i-1]=e[i];e[n-1]=0;var f=0;var tst1=0;var eps=1e-12;for(var l=0;l<n;l++){tst1=Math.max(tst1,Math.abs(d[l])+Math.abs(e[l]));var m=l;while(m<n){if(Math.abs(e[m])<=eps*tst1){break}m++}if(m>l){var iter=0;do{iter++;var g=d[l];var p=(d[l+1]-g)/(2*e[l]);var r=science.hypot(p,1);if(p<0)r=-r;d[l]=e[l]/(p+r);d[l+1]=e[l]*(p+r);var dl1=d[l+1];var h=g-d[l];for(var i=l+2;i<n;i++)d[i]-=h;f+=h;p=d[m];var c=1;var c2=c;var c3=c;var el1=e[l+1];var s=0;var s2=0;for(var i=m-1;i>=l;i--){c3=c2;c2=c;s2=s;g=c*e[i];h=c*p;r=science.hypot(p,e[i]);e[i+1]=s*r;s=e[i]/r;c=p/r;p=c*d[i]-s*g;d[i+1]=h+s*(c*g+s*d[i]);for(var k=0;k<n;k++){h=V[k][i+1];V[k][i+1]=s*V[k][i]+c*h;V[k][i]=c*V[k][i]-s*h}}p=-s*s2*c3*el1*e[l]/dl1;e[l]=s*p;d[l]=c*p}while(Math.abs(e[l])>eps*tst1)}d[l]=d[l]+f;e[l]=0}for(var i=0;i<n-1;i++){var k=i;var p=d[i];for(var j=i+1;j<n;j++){if(d[j]<p){k=j;p=d[j]}}if(k!=i){d[k]=d[i];d[i]=p;for(var j=0;j<n;j++){p=V[j][i];V[j][i]=V[j][k];V[j][k]=p}}}}function science_lin_decomposeOrthes(H,V){var n=H.length;var ort=[];var low=0;var high=n-1;for(var m=low+1;m<high;m++){var scale=0;for(var i=m;i<=high;i++)scale+=Math.abs(H[i][m-1]);if(scale!==0){var h=0;for(var i=high;i>=m;i--){ort[i]=H[i][m-1]/scale;h+=ort[i]*ort[i]}var g=Math.sqrt(h);if(ort[m]>0)g=-g;h=h-ort[m]*g;ort[m]=ort[m]-g;for(var j=m;j<n;j++){var f=0;for(var i=high;i>=m;i--)f+=ort[i]*H[i][j];f/=h;for(var i=m;i<=high;i++)H[i][j]-=f*ort[i]}for(var i=0;i<=high;i++){var f=0;for(var j=high;j>=m;j--)f+=ort[j]*H[i][j];f/=h;for(var j=m;j<=high;j++)H[i][j]-=f*ort[j]}ort[m]=scale*ort[m];H[m][m-1]=scale*g}}for(var i=0;i<n;i++){for(var j=0;j<n;j++)V[i][j]=i===j?1:0}for(var m=high-1;m>=low+1;m--){if(H[m][m-1]!==0){for(var i=m+1;i<=high;i++)ort[i]=H[i][m-1];for(var j=m;j<=high;j++){var g=0;for(var i=m;i<=high;i++)g+=ort[i]*V[i][j];g=g/ort[m]/H[m][m-1];for(var i=m;i<=high;i++)V[i][j]+=g*ort[i]}}}}function science_lin_decomposeHqr2(d,e,H,V){var nn=H.length,n=nn-1,low=0,high=nn-1,eps=1e-12,exshift=0,p=0,q=0,r=0,s=0,z=0,t,w,x,y;var norm=0;for(var i=0;i<nn;i++){if(i<low||i>high){d[i]=H[i][i];e[i]=0}for(var j=Math.max(i-1,0);j<nn;j++)norm+=Math.abs(H[i][j])}var iter=0;while(n>=low){var l=n;while(l>low){s=Math.abs(H[l-1][l-1])+Math.abs(H[l][l]);if(s===0)s=norm;if(Math.abs(H[l][l-1])<eps*s)break;l--}if(l===n){H[n][n]=H[n][n]+exshift;d[n]=H[n][n];e[n]=0;n--;iter=0}else if(l===n-1){w=H[n][n-1]*H[n-1][n];p=(H[n-1][n-1]-H[n][n])/2;q=p*p+w;z=Math.sqrt(Math.abs(q));H[n][n]=H[n][n]+exshift;H[n-1][n-1]=H[n-1][n-1]+exshift;x=H[n][n];if(q>=0){z=p+(p>=0?z:-z);d[n-1]=x+z;d[n]=d[n-1];if(z!==0)d[n]=x-w/z;e[n-1]=0;e[n]=0;x=H[n][n-1];s=Math.abs(x)+Math.abs(z);p=x/s;q=z/s;r=Math.sqrt(p*p+q*q);p/=r;q/=r;for(var j=n-1;j<nn;j++){z=H[n-1][j];H[n-1][j]=q*z+p*H[n][j];H[n][j]=q*H[n][j]-p*z}for(var i=0;i<=n;i++){z=H[i][n-1];H[i][n-1]=q*z+p*H[i][n];H[i][n]=q*H[i][n]-p*z}for(var i=low;i<=high;i++){z=V[i][n-1];V[i][n-1]=q*z+p*V[i][n];V[i][n]=q*V[i][n]-p*z}}else{d[n-1]=x+p;d[n]=x+p;e[n-1]=z;e[n]=-z}n=n-2;iter=0}else{x=H[n][n];y=0;w=0;if(l<n){y=H[n-1][n-1];w=H[n][n-1]*H[n-1][n]}if(iter==10){exshift+=x;for(var i=low;i<=n;i++){H[i][i]-=x}s=Math.abs(H[n][n-1])+Math.abs(H[n-1][n-2]);x=y=.75*s;w=-.4375*s*s}if(iter==30){s=(y-x)/2;s=s*s+w;if(s>0){s=Math.sqrt(s);if(y<x){s=-s}s=x-w/((y-x)/2+s);for(var i=low;i<=n;i++){H[i][i]-=s}exshift+=s;x=y=w=.964}}iter++;var m=n-2;while(m>=l){z=H[m][m];r=x-z;s=y-z;p=(r*s-w)/H[m+1][m]+H[m][m+1];q=H[m+1][m+1]-z-r-s;r=H[m+2][m+1];s=Math.abs(p)+Math.abs(q)+Math.abs(r);p=p/s;q=q/s;r=r/s;if(m==l)break;if(Math.abs(H[m][m-1])*(Math.abs(q)+Math.abs(r))<eps*(Math.abs(p)*(Math.abs(H[m-1][m-1])+Math.abs(z)+Math.abs(H[m+1][m+1])))){break}m--}for(var i=m+2;i<=n;i++){H[i][i-2]=0;if(i>m+2)H[i][i-3]=0}for(var k=m;k<=n-1;k++){var notlast=k!=n-1;if(k!=m){p=H[k][k-1];q=H[k+1][k-1];r=notlast?H[k+2][k-1]:0;x=Math.abs(p)+Math.abs(q)+Math.abs(r);if(x!=0){p/=x;q/=x;r/=x}}if(x==0)break;s=Math.sqrt(p*p+q*q+r*r);if(p<0){s=-s}if(s!=0){if(k!=m)H[k][k-1]=-s*x;else if(l!=m)H[k][k-1]=-H[k][k-1];p+=s;x=p/s;y=q/s;z=r/s;q/=p;r/=p;for(var j=k;j<nn;j++){p=H[k][j]+q*H[k+1][j];if(notlast){p=p+r*H[k+2][j];H[k+2][j]=H[k+2][j]-p*z}H[k][j]=H[k][j]-p*x;H[k+1][j]=H[k+1][j]-p*y}for(var i=0;i<=Math.min(n,k+3);i++){p=x*H[i][k]+y*H[i][k+1];if(notlast){p+=z*H[i][k+2];H[i][k+2]=H[i][k+2]-p*r}H[i][k]=H[i][k]-p;H[i][k+1]=H[i][k+1]-p*q}for(var i=low;i<=high;i++){p=x*V[i][k]+y*V[i][k+1];if(notlast){p=p+z*V[i][k+2];V[i][k+2]=V[i][k+2]-p*r}V[i][k]=V[i][k]-p;V[i][k+1]=V[i][k+1]-p*q}}}}}if(norm==0){return}for(n=nn-1;n>=0;n--){p=d[n];q=e[n];if(q==0){var l=n;H[n][n]=1;for(var i=n-1;i>=0;i--){w=H[i][i]-p;r=0;for(var j=l;j<=n;j++){r=r+H[i][j]*H[j][n]}if(e[i]<0){z=w;s=r}else{l=i;if(e[i]===0){H[i][n]=-r/(w!==0?w:eps*norm)}else{x=H[i][i+1];y=H[i+1][i];q=(d[i]-p)*(d[i]-p)+e[i]*e[i];t=(x*s-z*r)/q;H[i][n]=t;if(Math.abs(x)>Math.abs(z)){H[i+1][n]=(-r-w*t)/x}else{H[i+1][n]=(-s-y*t)/z}}t=Math.abs(H[i][n]);if(eps*t*t>1){for(var j=i;j<=n;j++)H[j][n]=H[j][n]/t}}}}else if(q<0){var l=n-1;if(Math.abs(H[n][n-1])>Math.abs(H[n-1][n])){H[n-1][n-1]=q/H[n][n-1];H[n-1][n]=-(H[n][n]-p)/H[n][n-1]}else{var zz=science_lin_decomposeCdiv(0,-H[n-1][n],H[n-1][n-1]-p,q);H[n-1][n-1]=zz[0];H[n-1][n]=zz[1]}H[n][n-1]=0;H[n][n]=1;for(var i=n-2;i>=0;i--){var ra=0,sa=0,vr,vi;for(var j=l;j<=n;j++){ra=ra+H[i][j]*H[j][n-1];sa=sa+H[i][j]*H[j][n]}w=H[i][i]-p;if(e[i]<0){z=w;r=ra;s=sa}else{l=i;if(e[i]==0){var zz=science_lin_decomposeCdiv(-ra,-sa,w,q);H[i][n-1]=zz[0];H[i][n]=zz[1]}else{x=H[i][i+1];y=H[i+1][i];vr=(d[i]-p)*(d[i]-p)+e[i]*e[i]-q*q;vi=(d[i]-p)*2*q;if(vr==0&vi==0){vr=eps*norm*(Math.abs(w)+Math.abs(q)+Math.abs(x)+Math.abs(y)+Math.abs(z))}var zz=science_lin_decomposeCdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi);H[i][n-1]=zz[0];H[i][n]=zz[1];if(Math.abs(x)>Math.abs(z)+Math.abs(q)){H[i+1][n-1]=(-ra-w*H[i][n-1]+q*H[i][n])/x;H[i+1][n]=(-sa-w*H[i][n]-q*H[i][n-1])/x}else{var zz=science_lin_decomposeCdiv(-r-y*H[i][n-1],-s-y*H[i][n],z,q);H[i+1][n-1]=zz[0];H[i+1][n]=zz[1]}}t=Math.max(Math.abs(H[i][n-1]),Math.abs(H[i][n]));if(eps*t*t>1){for(var j=i;j<=n;j++){H[j][n-1]=H[j][n-1]/t;H[j][n]=H[j][n]/t}}}}}}for(var i=0;i<nn;i++){if(i<low||i>high){for(var j=i;j<nn;j++)V[i][j]=H[i][j]}}for(var j=nn-1;j>=low;j--){for(var i=low;i<=high;i++){z=0;for(var k=low;k<=Math.min(j,high);k++)z+=V[i][k]*H[k][j];V[i][j]=z}}}function science_lin_decomposeCdiv(xr,xi,yr,yi){if(Math.abs(yr)>Math.abs(yi)){var r=yi/yr,d=yr+r*yi;return[(xr+r*xi)/d,(xi-r*xr)/d]}else{var r=yr/yi,d=yi+r*yr;return[(r*xr+xi)/d,(r*xi-xr)/d]}}science.lin.cross=function(a,b){return[a[1]*b[2]-a[2]*b[1],a[2]*b[0]-a[0]*b[2],a[0]*b[1]-a[1]*b[0]]};science.lin.dot=function(a,b){var s=0,i=-1,n=Math.min(a.length,b.length);while(++i<n)s+=a[i]*b[i];return s};science.lin.length=function(p){return Math.sqrt(science.lin.dot(p,p))};science.lin.normalize=function(p){var length=science.lin.length(p);return p.map(function(d){return d/length})};science.lin.determinant=function(matrix){var m=matrix[0].concat(matrix[1]).concat(matrix[2]).concat(matrix[3]);return m[12]*m[9]*m[6]*m[3]-m[8]*m[13]*m[6]*m[3]-m[12]*m[5]*m[10]*m[3]+m[4]*m[13]*m[10]*m[3]+m[8]*m[5]*m[14]*m[3]-m[4]*m[9]*m[14]*m[3]-m[12]*m[9]*m[2]*m[7]+m[8]*m[13]*m[2]*m[7]+m[12]*m[1]*m[10]*m[7]-m[0]*m[13]*m[10]*m[7]-m[8]*m[1]*m[14]*m[7]+m[0]*m[9]*m[14]*m[7]+m[12]*m[5]*m[2]*m[11]-m[4]*m[13]*m[2]*m[11]-m[12]*m[1]*m[6]*m[11]+m[0]*m[13]*m[6]*m[11]+m[4]*m[1]*m[14]*m[11]-m[0]*m[5]*m[14]*m[11]-m[8]*m[5]*m[2]*m[15]+m[4]*m[9]*m[2]*m[15]+m[8]*m[1]*m[6]*m[15]-m[0]*m[9]*m[6]*m[15]-m[4]*m[1]*m[10]*m[15]+m[0]*m[5]*m[10]*m[15]};science.lin.gaussjordan=function(m,eps){if(!eps)eps=1e-10;var h=m.length,w=m[0].length,y=-1,y2,x;while(++y<h){var maxrow=y;y2=y;while(++y2<h){if(Math.abs(m[y2][y])>Math.abs(m[maxrow][y]))maxrow=y2}var tmp=m[y];m[y]=m[maxrow];m[maxrow]=tmp;if(Math.abs(m[y][y])<=eps)return false;y2=y;while(++y2<h){var c=m[y2][y]/m[y][y];x=y-1;while(++x<w){m[y2][x]-=m[y][x]*c}}}y=h;while(--y>=0){var c=m[y][y];y2=-1;while(++y2<y){x=w;while(--x>=y){m[y2][x]-=m[y][x]*m[y2][y]/c}}m[y][y]/=c;x=h-1;while(++x<w){m[y][x]/=c}}return true};science.lin.inverse=function(m){var n=m.length,i=-1;if(n!==m[0].length)return;m=m.map(function(row,i){var identity=new Array(n),j=-1;while(++j<n)identity[j]=i===j?1:0;return row.concat(identity)});science.lin.gaussjordan(m);while(++i<n){m[i]=m[i].slice(n)}return m};science.lin.multiply=function(a,b){var m=a.length,n=b[0].length,p=b.length,i=-1,j,k;if(p!==a[0].length)throw{error:"columns(a) != rows(b); "+a[0].length+" != "+p};var ab=new Array(m);while(++i<m){ab[i]=new Array(n);j=-1;while(++j<n){var s=0;k=-1;while(++k<p)s+=a[i][k]*b[k][j];ab[i][j]=s}}return ab};science.lin.transpose=function(a){var m=a.length,n=a[0].length,i=-1,j,b=new Array(n);while(++i<n){b[i]=new Array(m);j=-1;while(++j<m)b[i][j]=a[j][i]}return b};science.lin.tridag=function(a,b,c,d,x,n){var i,m;for(i=1;i<n;i++){m=a[i]/b[i-1];b[i]-=m*c[i-1];d[i]-=m*d[i-1]}x[n-1]=d[n-1]/b[n-1];for(i=n-2;i>=0;i--){x[i]=(d[i]-c[i]*x[i+1])/b[i]}}})(this);(function(exports){science.stats={};science.stats.bandwidth={nrd0:function(x){var hi=Math.sqrt(science.stats.variance(x));if(!(lo=Math.min(hi,science.stats.iqr(x)/1.34)))(lo=hi)||(lo=Math.abs(x[1]))||(lo=1);return.9*lo*Math.pow(x.length,-.2)},nrd:function(x){var h=science.stats.iqr(x)/1.34;return 1.06*Math.min(Math.sqrt(science.stats.variance(x)),h)*Math.pow(x.length,-1/5)}};science.stats.distance={euclidean:function(a,b){var n=a.length,i=-1,s=0,x;while(++i<n){x=a[i]-b[i];s+=x*x}return Math.sqrt(s)},manhattan:function(a,b){var n=a.length,i=-1,s=0;while(++i<n)s+=Math.abs(a[i]-b[i]);return s},minkowski:function(p){return function(a,b){var n=a.length,i=-1,s=0;while(++i<n)s+=Math.pow(Math.abs(a[i]-b[i]),p);return Math.pow(s,1/p)}},chebyshev:function(a,b){var n=a.length,i=-1,max=0,x;while(++i<n){x=Math.abs(a[i]-b[i]);if(x>max)max=x}return max},hamming:function(a,b){var n=a.length,i=-1,d=0;while(++i<n)if(a[i]!==b[i])d++;return d},jaccard:function(a,b){var n=a.length,i=-1,s=0;while(++i<n)if(a[i]===b[i])s++;return s/n},braycurtis:function(a,b){var n=a.length,i=-1,s0=0,s1=0,ai,bi;while(++i<n){ai=a[i];bi=b[i];s0+=Math.abs(ai-bi);s1+=Math.abs(ai+bi)}return s0/s1}};science.stats.erf=function(x){var a1=.254829592,a2=-.284496736,a3=1.421413741,a4=-1.453152027,a5=1.061405429,p=.3275911;var sign=x<0?-1:1;if(x<0){sign=-1;x=-x}var t=1/(1+p*x);return sign*(1-((((a5*t+a4)*t+a3)*t+a2)*t+a1)*t*Math.exp(-x*x))};science.stats.phi=function(x){return.5*(1+science.stats.erf(x/Math.SQRT2))};science.stats.kernel={uniform:function(u){if(u<=1&&u>=-1)return.5;return 0},triangular:function(u){if(u<=1&&u>=-1)return 1-Math.abs(u);return 0},epanechnikov:function(u){if(u<=1&&u>=-1)return.75*(1-u*u);return 0},quartic:function(u){if(u<=1&&u>=-1){var tmp=1-u*u;return 15/16*tmp*tmp}return 0},triweight:function(u){if(u<=1&&u>=-1){var tmp=1-u*u;return 35/32*tmp*tmp*tmp}return 0},gaussian:function(u){return 1/Math.sqrt(2*Math.PI)*Math.exp(-.5*u*u)},cosine:function(u){if(u<=1&&u>=-1)return Math.PI/4*Math.cos(Math.PI/2*u);return 0}};science.stats.kde=function(){var kernel=science.stats.kernel.gaussian,sample=[],bandwidth=science.stats.bandwidth.nrd;function kde(points,i){var bw=bandwidth.call(this,sample);return points.map(function(x){var i=-1,y=0,n=sample.length;while(++i<n){y+=kernel((x-sample[i])/bw)}return[x,y/bw/n]})}kde.kernel=function(x){if(!arguments.length)return kernel;kernel=x;return kde};kde.sample=function(x){if(!arguments.length)return sample;sample=x;return kde};kde.bandwidth=function(x){if(!arguments.length)return bandwidth;bandwidth=science.functor(x);return kde};return kde};science.stats.kmeans=function(){var distance=science.stats.distance.euclidean,maxIterations=1e3,k=1;function kmeans(vectors){var n=vectors.length,assignments=[],clusterSizes=[],repeat=1,iterations=0,centroids=science_stats_kmeansRandom(k,vectors),newCentroids,i,j,x,d,min,best;while(repeat&&iterations<maxIterations){j=-1;while(++j<k){clusterSizes[j]=0}i=-1;while(++i<n){x=vectors[i];min=Infinity;j=-1;while(++j<k){d=distance.call(this,centroids[j],x);if(d<min){min=d;best=j}}clusterSizes[assignments[i]=best]++}newCentroids=[];i=-1;while(++i<n){x=assignments[i];d=newCentroids[x];if(d==null)newCentroids[x]=vectors[i].slice();else{j=-1;while(++j<d.length){d[j]+=vectors[i][j]}}}j=-1;while(++j<k){x=newCentroids[j];d=1/clusterSizes[j];i=-1;while(++i<x.length)x[i]*=d}repeat=0;j=-1;while(++j<k){if(!science_stats_kmeansCompare(newCentroids[j],centroids[j])){repeat=1;break}}centroids=newCentroids;iterations++}return{assignments:assignments,centroids:centroids}}kmeans.k=function(x){if(!arguments.length)return k;k=x;return kmeans};kmeans.distance=function(x){if(!arguments.length)return distance;distance=x;return kmeans};return kmeans};function science_stats_kmeansCompare(a,b){if(!a||!b||a.length!==b.length)return false;var n=a.length,i=-1;while(++i<n)if(a[i]!==b[i])return false;return true}function science_stats_kmeansRandom(k,vectors){var n=vectors.length;if(k>n)return null;var selected_vectors=[];var selected_indices=[];var tested_indices={};var tested=0;var selected=0;var i,vector,select;while(selected<k){if(tested===n)return null;var random_index=Math.floor(Math.random()*n);if(random_index in tested_indices)continue;tested_indices[random_index]=1;tested++;vector=vectors[random_index];select=true;for(i=0;i<selected;i++){if(science_stats_kmeansCompare(vector,selected_vectors[i])){select=false;break}}if(select){selected_vectors[selected]=vector;selected_indices[selected]=random_index;selected++}}return selected_vectors}science.stats.hcluster=function(){var distance=science.stats.distance.euclidean,linkage="single";function hcluster(vectors){var n=vectors.length,dMin=[],cSize=[],distMatrix=[],clusters=[],c1,c2,c1Cluster,c2Cluster,p,root,i,j;i=-1;while(++i<n){dMin[i]=0;distMatrix[i]=[];j=-1;while(++j<n){distMatrix[i][j]=i===j?Infinity:distance(vectors[i],vectors[j]);if(distMatrix[i][dMin[i]]>distMatrix[i][j])dMin[i]=j}}i=-1;while(++i<n){clusters[i]=[];clusters[i][0]={left:null,right:null,dist:0,centroid:vectors[i],size:1,depth:0};cSize[i]=1}for(p=0;p<n-1;p++){c1=0;for(i=0;i<n;i++){if(distMatrix[i][dMin[i]]<distMatrix[c1][dMin[c1]])c1=i}c2=dMin[c1];c1Cluster=clusters[c1][0];c2Cluster=clusters[c2][0];var newCluster={left:c1Cluster,right:c2Cluster,dist:distMatrix[c1][c2],centroid:calculateCentroid(c1Cluster.size,c1Cluster.centroid,c2Cluster.size,c2Cluster.centroid),size:c1Cluster.size+c2Cluster.size,depth:1+Math.max(c1Cluster.depth,c2Cluster.depth)};clusters[c1].splice(0,0,newCluster);cSize[c1]+=cSize[c2];for(j=0;j<n;j++){switch(linkage){case"single":if(distMatrix[c1][j]>distMatrix[c2][j])distMatrix[j][c1]=distMatrix[c1][j]=distMatrix[c2][j];break;case"complete":if(distMatrix[c1][j]<distMatrix[c2][j])distMatrix[j][c1]=distMatrix[c1][j]=distMatrix[c2][j];break;case"average":distMatrix[j][c1]=distMatrix[c1][j]=(cSize[c1]*distMatrix[c1][j]+cSize[c2]*distMatrix[c2][j])/(cSize[c1]+cSize[j]);break}}distMatrix[c1][c1]=Infinity;for(i=0;i<n;i++)distMatrix[i][c2]=distMatrix[c2][i]=Infinity;for(j=0;j<n;j++){if(dMin[j]==c2)dMin[j]=c1;if(distMatrix[c1][j]<distMatrix[c1][dMin[c1]])dMin[c1]=j}root=newCluster}return root}hcluster.distance=function(x){if(!arguments.length)return distance;distance=x;return hcluster};return hcluster};function calculateCentroid(c1Size,c1Centroid,c2Size,c2Centroid){var newCentroid=[],newSize=c1Size+c2Size,n=c1Centroid.length,i=-1;while(++i<n){newCentroid[i]=(c1Size*c1Centroid[i]+c2Size*c2Centroid[i])/newSize}return newCentroid}science.stats.iqr=function(x){var quartiles=science.stats.quantiles(x,[.25,.75]);return quartiles[1]-quartiles[0]};science.stats.loess=function(){var bandwidth=.3,robustnessIters=2,accuracy=1e-12;function smooth(xval,yval,weights){var n=xval.length,i;if(n!==yval.length)throw{error:"Mismatched array lengths"};if(n==0)throw{error:"At least one point required."};if(arguments.length<3){weights=[];i=-1;while(++i<n)weights[i]=1}science_stats_loessFiniteReal(xval);science_stats_loessFiniteReal(yval);science_stats_loessFiniteReal(weights);science_stats_loessStrictlyIncreasing(xval);if(n==1)return[yval[0]];if(n==2)return[yval[0],yval[1]];var bandwidthInPoints=Math.floor(bandwidth*n);if(bandwidthInPoints<2)throw{error:"Bandwidth too small."};var res=[],residuals=[],robustnessWeights=[];i=-1;while(++i<n){res[i]=0;residuals[i]=0;robustnessWeights[i]=1}var iter=-1;while(++iter<=robustnessIters){var bandwidthInterval=[0,bandwidthInPoints-1];var x;i=-1;while(++i<n){x=xval[i];if(i>0){science_stats_loessUpdateBandwidthInterval(xval,weights,i,bandwidthInterval)}var ileft=bandwidthInterval[0],iright=bandwidthInterval[1];var edge=xval[i]-xval[ileft]>xval[iright]-xval[i]?ileft:iright;var sumWeights=0,sumX=0,sumXSquared=0,sumY=0,sumXY=0,denom=Math.abs(1/(xval[edge]-x));for(var k=ileft;k<=iright;++k){var xk=xval[k],yk=yval[k],dist=k<i?x-xk:xk-x,w=science_stats_loessTricube(dist*denom)*robustnessWeights[k]*weights[k],xkw=xk*w;sumWeights+=w;sumX+=xkw;sumXSquared+=xk*xkw;sumY+=yk*w;sumXY+=yk*xkw}var meanX=sumX/sumWeights,meanY=sumY/sumWeights,meanXY=sumXY/sumWeights,meanXSquared=sumXSquared/sumWeights;var beta=Math.sqrt(Math.abs(meanXSquared-meanX*meanX))<accuracy?0:(meanXY-meanX*meanY)/(meanXSquared-meanX*meanX);var alpha=meanY-beta*meanX;res[i]=beta*x+alpha;residuals[i]=Math.abs(yval[i]-res[i])}if(iter===robustnessIters){break}var medianResidual=science.stats.median(residuals);if(Math.abs(medianResidual)<accuracy)break;var arg,w;i=-1;while(++i<n){arg=residuals[i]/(6*medianResidual);robustnessWeights[i]=arg>=1?0:(w=1-arg*arg)*w}}return res}smooth.bandwidth=function(x){if(!arguments.length)return x;bandwidth=x;return smooth};smooth.robustnessIterations=function(x){if(!arguments.length)return x;robustnessIters=x;return smooth};smooth.accuracy=function(x){if(!arguments.length)return x;accuracy=x;return smooth};return smooth};function science_stats_loessFiniteReal(values){var n=values.length,i=-1;while(++i<n)if(!isFinite(values[i]))return false;return true}function science_stats_loessStrictlyIncreasing(xval){var n=xval.length,i=0;while(++i<n)if(xval[i-1]>=xval[i])return false;return true}function science_stats_loessTricube(x){return(x=1-x*x*x)*x*x}function science_stats_loessUpdateBandwidthInterval(xval,weights,i,bandwidthInterval){var left=bandwidthInterval[0],right=bandwidthInterval[1];var nextRight=science_stats_loessNextNonzero(weights,right);if(nextRight<xval.length&&xval[nextRight]-xval[i]<xval[i]-xval[left]){var nextLeft=science_stats_loessNextNonzero(weights,left);bandwidthInterval[0]=nextLeft;bandwidthInterval[1]=nextRight}}function science_stats_loessNextNonzero(weights,i){var j=i+1;while(j<weights.length&&weights[j]===0)j++;return j}science.stats.mean=function(x){var n=x.length;if(n===0)return NaN;var m=0,i=-1;while(++i<n)m+=(x[i]-m)/(i+1);return m};science.stats.median=function(x){return science.stats.quantiles(x,[.5])[0]};science.stats.mode=function(x){var counts={},mode=[],max=0,n=x.length,i=-1,d,k;while(++i<n){k=counts.hasOwnProperty(d=x[i])?++counts[d]:counts[d]=1;if(k===max)mode.push(d);else if(k>max){max=k;mode=[d]}}if(mode.length===1)return mode[0]};science.stats.quantiles=function(d,quantiles){d=d.slice().sort(science.ascending);var n_1=d.length-1;return quantiles.map(function(q){if(q===0)return d[0];else if(q===1)return d[n_1];var index=1+q*n_1,lo=Math.floor(index),h=index-lo,a=d[lo-1];return h===0?a:a+h*(d[lo]-a)})};science.stats.variance=function(x){var n=x.length;if(n<1)return NaN;if(n===1)return 0;var mean=science.stats.mean(x),i=-1,s=0;while(++i<n){var v=x[i]-mean;s+=v*v}return s/(n-1)};science.stats.distribution={};science.stats.distribution.gaussian=function(){var random=Math.random,mean=0,sigma=1,variance=1;function gaussian(){var x1,x2,rad,y1;do{x1=2*random()-1;x2=2*random()-1;rad=x1*x1+x2*x2}while(rad>=1||rad===0);return mean+sigma*x1*Math.sqrt(-2*Math.log(rad)/rad)}gaussian.pdf=function(x){x=(x-mean)/sigma;return science_stats_distribution_gaussianConstant*Math.exp(-.5*x*x)/sigma};gaussian.cdf=function(x){x=(x-mean)/sigma;return.5*(1+science.stats.erf(x/Math.SQRT2))};gaussian.mean=function(x){if(!arguments.length)return mean;mean=+x;return gaussian};gaussian.variance=function(x){if(!arguments.length)return variance;sigma=Math.sqrt(variance=+x);return gaussian};gaussian.random=function(x){if(!arguments.length)return random;random=x;return gaussian};return gaussian};science_stats_distribution_gaussianConstant=1/Math.sqrt(2*Math.PI)})(this)})(this);
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment