Skip to content

Instantly share code, notes, and snippets.

@amergin
Last active August 29, 2015 14:26
Show Gist options
  • Save amergin/0e88075d9d94141c3094 to your computer and use it in GitHub Desktop.
Save amergin/0e88075d9d94141c3094 to your computer and use it in GitHub Desktop.
Math.js regression performance problem, code
_xMatrixTransp size = 2,4084
_xMatrix size = 4084
_normalTargetData size = 4084
function numericJs() {
var xMatrix = numeric.transpose(xMatrixTransp);
// see https://en.wikipedia.org/wiki/Ordinary_least_squares#Estimation
// Compute beta = (X^T X)^{-1} X^T y
var dotProduct = numeric.dot(xMatrixTransp, xMatrix);
var inverse = numeric.inv(dotProduct);
dotProduct = null;
var multi2 = numeric.dot(inverse, xMatrixTransp);
var betas = numeric.dot(multi2, normalTargetData);
multi2 = null;
// console.log("numericjs betas", betas);
var n = _.size(xMatrix),
k = _.size(xMatrix[0]),
beta = betas[1];
// console.log("n, k", n, k);
// get confidence interval and p-value
var ciAndPvalue = regressionUtils.getCIAndPvalueNumericjs(inverse, xMatrix, xMatrixTransp, [normalTargetData], n, k, beta);
return {
result: { success: true },
ci: ciAndPvalue.ci,
pvalue: ciAndPvalue.pvalue,
betas: betas
};
}
function mathJs() {
// see https://en.wikipedia.org/wiki/Ordinary_least_squares#Estimation
// beta = (X^T X)^{-1} X^T y
var _xMatrixTransp = math.matrix(xMatrixTransp);
var _xMatrix = math.transpose(xMatrixTransp);
var _normalTargetData = math.matrix(normalTargetData);
console.log("_xMatrixTransp size = ", String(_xMatrixTransp.size()));
console.log("_xMatrix size = ", _.size(_xMatrix));
console.log("_normalTargetData size = ", String(_normalTargetData.size()));
// Compute beta = (X^T X)^{-1} X^T y
var dotProduct = math.multiply(_xMatrixTransp, _xMatrix);
var inverse = math.inv(dotProduct);
var multi2 = math.multiply(inverse, _xMatrixTransp);
var betas = math.multiply(multi2, normalTargetData);
// console.log("mathjs betas=", String(betas));
var n = math.size(_xMatrix)[0];
var k = math.size(_xMatrix)[1];
var beta = betas.subset(math.index(1));
var _yMatrixTransp = math.matrix([normalTargetData]);
// get confidence interval and p-value
var ciAndPvalue = regressionUtils.getCIAndPvalueMathjs(inverse, _xMatrix, _xMatrixTransp, _yMatrixTransp, n, k, beta);
return {
result: { success: true },
ci: ciAndPvalue.ci,
pvalue: ciAndPvalue.pvalue,
betas: betas.valueOf()
};
}
root.getCIAndPvalueMathjs = function(dotInverse, xMatrix, xMatrixTransposed, yMatrixTransposed, n, k, beta) {
var VARIABLE_INDEX = 1;
// H = X (X^T X)^{-1} X^T
var hMatrix = math.chain(xMatrix)
.multiply(dotInverse)
.multiply(xMatrixTransposed)
.done();
xMatrix = null;
var hMatrixSize = math.size(hMatrix).subset(math.index(0));
var identity = math.eye(hMatrixSize, 'sparse');
var subtracted = math.subtract(identity, hMatrix);
var yMatrix = math.transpose(yMatrixTransposed);
var degrees = n-(k+1);
var sigma = math.chain(yMatrixTransposed)
.multiply(subtracted)
.multiply(yMatrix)
.divide(degrees)
.done();
yMatrix = null;
subtracted = null;
identity = null;
yMatrixTransposed = null;
var cMatrix = math.multiply(sigma, dotInverse);
dotInverse = null;
var alpha = 0.05 / k;
var sqrt = Math.sqrt( cMatrix.subset(math.index(VARIABLE_INDEX, VARIABLE_INDEX)) );
cMatrix = null;
var ci = statDist.tdistr(degrees, alpha/2) * sqrt;
// See http://reliawiki.org/index.php/Multiple_Linear_Regression_Analysis
// -> p value = 2 * (1-P(T <= |t0|)
var t = beta / sqrt;
sqrt = null;
var pvalue = statDist.tprob(degrees, t);
return {
ci: [ beta - ci, beta + ci ],
pvalue: pvalue
};
};
root.getCIAndPvalueNumericjs = function(dotInverse, xMatrix, xMatrixTransposed, yMatrixTransposed, n, k, beta) {
var VARIABLE_INDEX = 1;
var hMatrix = numeric.dot( numeric.dot(xMatrix, dotInverse), xMatrixTransposed );
var _identity = numeric.identity( _.size(hMatrix) );
var _subtracted = numeric.sub(_identity, hMatrix);
var yMatrix = numeric.transpose(yMatrixTransposed);
var sigma = numeric.dot( numeric.dot( yMatrixTransposed, _subtracted ), yMatrix )[0][0] / (n-(k+1));
var cMatrix = numeric.mul(sigma, dotInverse);
var alpha = 0.05 / k;
var _sqrt = Math.sqrt( cMatrix[VARIABLE_INDEX][VARIABLE_INDEX] ); //numeric.getDiag(cMatrix)[1] );
var degrees = n-(k+1);
var ci = statDist.tdistr(degrees, alpha/2) * _sqrt;
// See http://reliawiki.org/index.php/Multiple_Linear_Regression_Analysis
// -> p value = 2 * (1-P(T <= |t0|)
var t = beta / _sqrt;
var pvalue = statDist.tprob(degrees, t);
return {
ci: [ beta - ci, beta + ci ],
pvalue: pvalue
};
};
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment