Instantly share code, notes, and snippets.

# amergin/output.txt

Last active August 29, 2015 14:26
Show Gist options
• Save amergin/0e88075d9d94141c3094 to your computer and use it in GitHub Desktop.
Math.js regression performance problem, code
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 _xMatrixTransp size = 2,4084 _xMatrix size = 4084 _normalTargetData size = 4084
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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 }; };