this.j$ = this.jStat = (function(Math, undefined) { // For quick reference. var concat = Array.prototype.concat; var slice = Array.prototype.slice; var toString = Object.prototype.toString; // Calculate correction for IEEE error // TODO: This calculation can be improved. function calcRdx(n, m) { var val = n > m ? n : m; return Math.pow(10, 17 - ~~(Math.log(((val > 0) ? val : -val)) * Math.LOG10E)); } var isArray = Array.isArray || function isArray(arg) { return toString.call(arg) === '[object Array]'; }; function isFunction(arg) { return toString.call(arg) === '[object Function]'; } function isNumber(arg) { return typeof arg === 'number' && arg === arg; } // Converts the jStat matrix to vector. function toVector(arr) { return concat.apply([], arr); } // The one and only jStat constructor. function jStat() { return new jStat._init(arguments); } // TODO: Remove after all references in src files have been removed. jStat.fn = jStat.prototype; // By separating the initializer from the constructor it's easier to handle // always returning a new instance whether "new" was used or not. jStat._init = function _init(args) { var i; // If first argument is an array, must be vector or matrix. if (isArray(args[0])) { // Check if matrix. if (isArray(args[0][0])) { // See if a mapping function was also passed. if (isFunction(args[1])) args[0] = jStat.map(args[0], args[1]); // Iterate over each is faster than this.push.apply(this, args[0]. for (i = 0; i < args[0].length; i++) this[i] = args[0][i]; this.length = args[0].length; // Otherwise must be a vector. } else { this[0] = isFunction(args[1]) ? jStat.map(args[0], args[1]) : args[0]; this.length = 1; } // If first argument is number, assume creation of sequence. } else if (isNumber(args[0])) { this[0] = jStat.seq.apply(null, args); this.length = 1; // Handle case when jStat object is passed to jStat. } else if (args[0] instanceof jStat) { // Duplicate the object and pass it back. return jStat(args[0].toArray()); // Unexpected argument value, return empty jStat object. // TODO: This is strange behavior. Shouldn't this throw or some such to let // the user know they had bad arguments? } else { this[0] = []; this.length = 1; } return this; }; jStat._init.prototype = jStat.prototype; jStat._init.constructor = jStat; // Utility functions. // TODO: for internal use only? jStat.utils = { calcRdx: calcRdx, isArray: isArray, isFunction: isFunction, isNumber: isNumber, toVector: toVector }; // Easily extend the jStat object. // TODO: is this seriously necessary? jStat.extend = function extend(obj) { var i, j; if (arguments.length === 1) { for (j in obj) jStat[j] = obj[j]; return this; } for (i = 1; i < arguments.length; i++) { for (j in arguments[i]) obj[j] = arguments[i][j]; } return obj; }; // Returns the number of rows in the matrix. jStat.rows = function rows(arr) { return arr.length || 1; }; // Returns the number of columns in the matrix. jStat.cols = function cols(arr) { return arr[0].length || 1; }; // Returns the dimensions of the object { rows: i, cols: j } jStat.dimensions = function dimensions(arr) { return { rows: jStat.rows(arr), cols: jStat.cols(arr) }; }; // Returns a specified row as a vector jStat.row = function row(arr, index) { return arr[index]; }; // Returns the specified column as a vector jStat.col = function cols(arr, index) { var column = new Array(arr.length); for (var i = 0; i < arr.length; i++) column[i] = [arr[i][index]]; return column; }; // Returns the diagonal of the matrix jStat.diag = function diag(arr) { var nrow = jStat.rows(arr); var res = new Array(nrow); for (var row = 0; row < nrow; row++) res[row] = [arr[row][row]]; return res; }; // Returns the anti-diagonal of the matrix jStat.antidiag = function antidiag(arr) { var nrow = jStat.rows(arr) - 1; var res = new Array(nrow); for (var i = 0; nrow >= 0; nrow--, i++) res[i] = [arr[i][nrow]]; return res; }; // Transpose a matrix or array. jStat.transpose = function transpose(arr) { var obj = []; var objArr, rows, cols, j, i; // Make sure arr is in matrix format. if (!isArray(arr[0])) arr = [arr]; rows = arr.length; cols = arr[0].length; for (i = 0; i < cols; i++) { objArr = new Array(rows); for (j = 0; j < rows; j++) objArr[j] = arr[j][i]; obj.push(objArr); } // If obj is vector, return only single array. return obj.length === 1 ? obj[0] : obj; }; // Map a function to an array or array of arrays. // "toAlter" is an internal variable. jStat.map = function map(arr, func, toAlter) { var row, nrow, ncol, res, col; if (!isArray(arr[0])) arr = [arr]; nrow = arr.length; ncol = arr[0].length; res = toAlter ? arr : new Array(nrow); for (row = 0; row < nrow; row++) { // if the row doesn't exist, create it if (!res[row]) res[row] = new Array(ncol); for (col = 0; col < ncol; col++) res[row][col] = func(arr[row][col], row, col); } return res.length === 1 ? res[0] : res; }; // Cumulatively combine the elements of an array or array of arrays using a function. jStat.cumreduce = function cumreduce(arr, func, toAlter) { var row, nrow, ncol, res, col; if (!isArray(arr[0])) arr = [arr]; nrow = arr.length; ncol = arr[0].length; res = toAlter ? arr : new Array(nrow); for (row = 0; row < nrow; row++) { // if the row doesn't exist, create it if (!res[row]) res[row] = new Array(ncol); if (ncol > 0) res[row][0] = arr[row][0]; for (col = 1; col < ncol; col++) res[row][col] = func(res[row][col-1], arr[row][col]); } return res.length === 1 ? res[0] : res; }; // Destructively alter an array. jStat.alter = function alter(arr, func) { return jStat.map(arr, func, true); }; // Generate a rows x cols matrix according to the supplied function. jStat.create = function create(rows, cols, func) { var res = new Array(rows); var i, j; if (isFunction(cols)) { func = cols; cols = rows; } for (i = 0; i < rows; i++) { res[i] = new Array(cols); for (j = 0; j < cols; j++) res[i][j] = func(i, j); } return res; }; function retZero() { return 0; } // Generate a rows x cols matrix of zeros. jStat.zeros = function zeros(rows, cols) { if (!isNumber(cols)) cols = rows; return jStat.create(rows, cols, retZero); }; function retOne() { return 1; } // Generate a rows x cols matrix of ones. jStat.ones = function ones(rows, cols) { if (!isNumber(cols)) cols = rows; return jStat.create(rows, cols, retOne); }; // Generate a rows x cols matrix of uniformly random numbers. jStat.rand = function rand(rows, cols) { if (!isNumber(cols)) cols = rows; return jStat.create(rows, cols, Math.random); }; function retIdent(i, j) { return i === j ? 1 : 0; } // Generate an identity matrix of size row x cols. jStat.identity = function identity(rows, cols) { if (!isNumber(cols)) cols = rows; return jStat.create(rows, cols, retIdent); }; // Tests whether a matrix is symmetric jStat.symmetric = function symmetric(arr) { var issymmetric = true; var size = arr.length; var row, col; if (arr.length !== arr[0].length) return false; for (row = 0; row < size; row++) { for (col = 0; col < size; col++) if (arr[col][row] !== arr[row][col]) return false; } return true; }; // Set all values to zero. jStat.clear = function clear(arr) { return jStat.alter(arr, retZero); }; // Generate sequence. jStat.seq = function seq(min, max, length, func) { if (!isFunction(func)) func = false; var arr = []; var hival = calcRdx(min, max); var step = (max * hival - min * hival) / ((length - 1) * hival); var current = min; var cnt; // Current is assigned using a technique to compensate for IEEE error. // TODO: Needs better implementation. for (cnt = 0; current <= max; cnt++, current = (min * hival + step * hival * cnt) / hival) { arr.push((func ? func(current, cnt) : current)); } return arr; }; // TODO: Go over this entire implementation. Seems a tragic waste of resources // doing all this work. Instead, and while ugly, use new Function() to generate // a custom function for each static method. // Quick reference. var jProto = jStat.prototype; // Default length. jProto.length = 0; // For internal use only. // TODO: Check if they're actually used, and if they are then rename them // to _* jProto.push = Array.prototype.push; jProto.sort = Array.prototype.sort; jProto.splice = Array.prototype.splice; jProto.slice = Array.prototype.slice; // Return a clean array. jProto.toArray = function toArray() { return this.length > 1 ? slice.call(this) : slice.call(this)[0]; }; // Map a function to a matrix or vector. jProto.map = function map(func, toAlter) { return jStat(jStat.map(this, func, toAlter)); }; // Cumulatively combine the elements of a matrix or vector using a function. jProto.cumreduce = function cumreduce(func, toAlter) { return jStat(jStat.cumreduce(this, func, toAlter)); }; // Destructively alter an array. jProto.alter = function alter(func) { jStat.alter(this, func); return this; }; // Extend prototype with methods that have no argument. (function(funcs) { for (var i = 0; i < funcs.length; i++) (function(passfunc) { jProto[passfunc] = function(func) { var self = this, results; // Check for callback. if (func) { setTimeout(function() { func.call(self, jProto[passfunc].call(self)); }); return this; } results = jStat[passfunc](this); return isArray(results) ? jStat(results) : results; }; })(funcs[i]); })('transpose clear symmetric rows cols dimensions diag antidiag'.split(' ')); // Extend prototype with methods that have one argument. (function(funcs) { for (var i = 0; i < funcs.length; i++) (function(passfunc) { jProto[passfunc] = function(index, func) { var self = this; // check for callback if (func) { setTimeout(function() { func.call(self, jProto[passfunc].call(self, index)); }); return this; } return jStat(jStat[passfunc](this, index)); }; })(funcs[i]); })('row col'.split(' ')); // Extend prototype with simple shortcut methods. (function(funcs) { for (var i = 0; i < funcs.length; i++) (function(passfunc) { jProto[passfunc] = new Function( 'return jStat(jStat.' + passfunc + '.apply(null, arguments));'); })(funcs[i]); })('create zeros ones rand identity'.split(' ')); // Exposing jStat. return jStat; }(Math)); (function(jStat, Math) { var isFunction = jStat.utils.isFunction; // Ascending functions for sort function ascNum(a, b) { return a - b; } function clip(arg, min, max) { return Math.max(min, Math.min(arg, max)); } // sum of an array jStat.sum = function sum(arr) { var sum = 0; var i = arr.length; var tmp; while (--i >= 0) sum += arr[i]; return sum; }; // sum squared jStat.sumsqrd = function sumsqrd(arr) { var sum = 0; var i = arr.length; while (--i >= 0) sum += arr[i] * arr[i]; return sum; }; // sum of squared errors of prediction (SSE) jStat.sumsqerr = function sumsqerr(arr) { var mean = jStat.mean(arr); var sum = 0; var i = arr.length; var tmp; while (--i >= 0) { tmp = arr[i] - mean; sum += tmp * tmp; } return sum; }; // product of an array jStat.product = function product(arr) { var prod = 1; var i = arr.length; while (--i >= 0) prod *= arr[i]; return prod; }; // minimum value of an array jStat.min = function min(arr) { var low = arr[0]; var i = 0; while (++i < arr.length) if (arr[i] < low) low = arr[i]; return low; }; // maximum value of an array jStat.max = function max(arr) { var high = arr[0]; var i = 0; while (++i < arr.length) if (arr[i] > high) high = arr[i]; return high; }; // mean value of an array jStat.mean = function mean(arr) { return jStat.sum(arr) / arr.length; }; // mean squared error (MSE) jStat.meansqerr = function meansqerr(arr) { return jStat.sumsqerr(arr) / arr.length; }; // geometric mean of an array jStat.geomean = function geomean(arr) { return Math.pow(jStat.product(arr), 1 / arr.length); }; // median of an array jStat.median = function median(arr) { var arrlen = arr.length; var _arr = arr.slice().sort(ascNum); // check if array is even or odd, then return the appropriate return !(arrlen & 1) ? (_arr[(arrlen / 2) - 1 ] + _arr[(arrlen / 2)]) / 2 : _arr[(arrlen / 2) | 0 ]; }; // cumulative sum of an array jStat.cumsum = function cumsum(arr) { return jStat.cumreduce(arr, function (a, b) { return a + b; }); }; // cumulative product of an array jStat.cumprod = function cumprod(arr) { return jStat.cumreduce(arr, function (a, b) { return a * b; }); }; // successive differences of a sequence jStat.diff = function diff(arr) { var diffs = []; var arrLen = arr.length; var i; for (i = 1; i < arrLen; i++) diffs.push(arr[i] - arr[i - 1]); return diffs; }; // mode of an array // if there are multiple modes of an array, return all of them // is this the appropriate way of handling it? jStat.mode = function mode(arr) { var arrLen = arr.length; var _arr = arr.slice().sort(ascNum); var count = 1; var maxCount = 0; var numMaxCount = 0; var mode_arr = []; var i; for (i = 0; i < arrLen; i++) { if (_arr[i] === _arr[i + 1]) { count++; } else { if (count > maxCount) { mode_arr = [_arr[i]]; maxCount = count; numMaxCount = 0; } // are there multiple max counts else if (count === maxCount) { mode_arr.push(_arr[i]); numMaxCount++; } // resetting count for new value in array count = 1; } } return numMaxCount === 0 ? mode_arr[0] : mode_arr; }; // range of an array jStat.range = function range(arr) { return jStat.max(arr) - jStat.min(arr); }; // variance of an array // flag indicates population vs sample jStat.variance = function variance(arr, flag) { return jStat.sumsqerr(arr) / (arr.length - (flag ? 1 : 0)); }; // standard deviation of an array // flag indicates population vs sample jStat.stdev = function stdev(arr, flag) { return Math.sqrt(jStat.variance(arr, flag)); }; // mean deviation (mean absolute deviation) of an array jStat.meandev = function meandev(arr) { var devSum = 0; var mean = jStat.mean(arr); var i; for (i = arr.length - 1; i >= 0; i--) devSum += Math.abs(arr[i] - mean); return devSum / arr.length; }; // median deviation (median absolute deviation) of an array jStat.meddev = function meddev(arr) { var devSum = 0; var median = jStat.median(arr); var i; for (i = arr.length - 1; i >= 0; i--) devSum += Math.abs(arr[i] - median); return devSum / arr.length; }; // coefficient of variation jStat.coeffvar = function coeffvar(arr) { return jStat.stdev(arr) / jStat.mean(arr); }; // quartiles of an array jStat.quartiles = function quartiles(arr) { var arrlen = arr.length; var _arr = arr.slice().sort(ascNum); return [ _arr[ Math.round((arrlen) / 4) - 1 ], _arr[ Math.round((arrlen) / 2) - 1 ], _arr[ Math.round((arrlen) * 3 / 4) - 1 ] ]; }; // Arbitary quantiles of an array. Direct port of the scipy.stats // implementation by Pierre GF Gerard-Marchant. jStat.quantiles = function quantiles(arr, quantilesArray, alphap, betap) { var sortedArray = arr.slice().sort(ascNum); var quantileVals = [quantilesArray.length]; var n = arr.length; var i, p, m, aleph, k, gamma; if (typeof alphap === 'undefined') alphap = 3 / 8; if (typeof betap === 'undefined') betap = 3 / 8; for (i = 0; i < quantilesArray.length; i++) { p = quantilesArray[i]; m = alphap + p * (1 - alphap - betap); aleph = n * p + m; k = Math.floor(clip(aleph, 1, n - 1)); gamma = clip(aleph - k, 0, 1); quantileVals[i] = (1 - gamma) * sortedArray[k - 1] + gamma * sortedArray[k]; } return quantileVals; }; // The percentile rank of score in a given array. Returns the percentage // of all values in the input array that are less than (kind='strict') or // less or equal than (kind='weak') score. Default is weak. jStat.percentileOfScore = function percentileOfScore(arr, score, kind) { var counter = 0; var len = arr.length; var strict = false; var value, i; if (kind === 'strict') strict = true; for (i = 0; i < len; i++) { value = arr[i]; if ((strict && value < score) || (!strict && value <= score)) { counter++; } } return counter / len; }; // Histogram (bin count) data jStat.histogram = function histogram(arr, bins) { var first = jStat.min(arr); var binCnt = bins || 4; var binWidth = (jStat.max(arr) - first) / binCnt; var len = arr.length; var bins = []; var i; for (i = 0; i < binCnt; i++) bins[i] = 0; for (i = 0; i < len; i++) bins[Math.min(Math.floor(((arr[i] - first) / binWidth)), binCnt - 1)] += 1; return bins; }; // covariance of two arrays jStat.covariance = function covariance(arr1, arr2) { var u = jStat.mean(arr1); var v = jStat.mean(arr2); var arr1Len = arr1.length; var sq_dev = new Array(arr1Len); var i; for (i = 0; i < arr1Len; i++) sq_dev[i] = (arr1[i] - u) * (arr2[i] - v); return jStat.sum(sq_dev) / (arr1Len - 1); }; // (pearson's) population correlation coefficient, rho jStat.corrcoeff = function corrcoeff(arr1, arr2) { return jStat.covariance(arr1, arr2) / jStat.stdev(arr1, 1) / jStat.stdev(arr2, 1); }; // statistical standardized moments (general form of skew/kurt) jStat.stanMoment = function stanMoment(arr, n) { var mu = jStat.mean(arr); var sigma = jStat.stdev(arr); var len = arr.length; var skewSum = 0; for (i = 0; i < len; i++) skewSum += Math.pow((arr[i] - mu) / sigma, n); return skewSum / arr.length; }; // (pearson's) moment coefficient of skewness jStat.skewness = function skewness(arr) { return jStat.stanMoment(arr, 3); }; // (pearson's) (excess) kurtosis jStat.kurtosis = function kurtosis(arr) { return jStat.stanMoment(arr, 4) - 3; }; var jProto = jStat.prototype; // Extend jProto with method for calculating cumulative sums and products. // This differs from the similar extension below as cumsum and cumprod should // not be run again in the case fullbool === true. // If a matrix is passed, automatically assume operation should be done on the // columns. (function(funcs) { for (var i = 0; i < funcs.length; i++) (function(passfunc) { // If a matrix is passed, automatically assume operation should be done on // the columns. jProto[passfunc] = function(fullbool, func) { var arr = []; var i = 0; var tmpthis = this; // Assignment reassignation depending on how parameters were passed in. if (isFunction(fullbool)) { func = fullbool; fullbool = false; } // Check if a callback was passed with the function. if (func) { setTimeout(function() { func.call(tmpthis, jProto[passfunc].call(tmpthis, fullbool)); }); return this; } // Check if matrix and run calculations. if (this.length > 1) { tmpthis = fullbool === true ? this : this.transpose(); for (; i < tmpthis.length; i++) arr[i] = jStat[passfunc](tmpthis[i]); return arr; } // Pass fullbool if only vector, not a matrix. for variance and stdev. return jStat[passfunc](this[0], fullbool); }; })(funcs[i]); })(('cumsum cumprod').split(' ')); // Extend jProto with methods which don't require arguments and work on columns. (function(funcs) { for (var i = 0; i < funcs.length; i++) (function(passfunc) { // If a matrix is passed, automatically assume operation should be done on // the columns. jProto[passfunc] = function(fullbool, func) { var arr = []; var i = 0; var tmpthis = this; // Assignment reassignation depending on how parameters were passed in. if (isFunction(fullbool)) { func = fullbool; fullbool = false; } // Check if a callback was passed with the function. if (func) { setTimeout(function() { func.call(tmpthis, jProto[passfunc].call(tmpthis, fullbool)); }); return this; } // Check if matrix and run calculations. if (this.length > 1) { tmpthis = fullbool === true ? this : this.transpose(); for (; i < tmpthis.length; i++) arr[i] = jStat[passfunc](tmpthis[i]); return fullbool === true ? jStat[passfunc](jStat.utils.toVector(arr)) : arr; } // Pass fullbool if only vector, not a matrix. for variance and stdev. return jStat[passfunc](this[0], fullbool); }; })(funcs[i]); })(('sum sumsqrd sumsqerr product min max mean meansqerr geomean median diff ' + 'mode range variance stdev meandev meddev coeffvar quartiles histogram ' + 'skewness kurtosis').split(' ')); // Extend jProto with functions that take arguments. Operations on matrices are // done on columns. (function(funcs) { for (var i = 0; i < funcs.length; i++) (function(passfunc) { jProto[passfunc] = function() { var arr = []; var i = 0; var tmpthis = this; var args = Array.prototype.slice.call(arguments); // If the last argument is a function, we assume it's a callback; we // strip the callback out and call the function again. if (isFunction(args[args.length - 1])) { var callbackFunction = args[args.length - 1]; var argsToPass = args.slice(0, args.length - 1); setTimeout(function() { callbackFunction.call(tmpthis, jProto[passfunc].apply(tmpthis, argsToPass)); }); return this; // Otherwise we curry the function args and call normally. } else { var callbackFunction = undefined; var curriedFunction = function curriedFunction(vector) { return jStat[passfunc].apply(tmpthis, [vector].concat(args)); } } // If this is a matrix, run column-by-column. if (this.length > 1) { tmpthis = tmpthis.transpose(); for (; i < tmpthis.length; i++) arr[i] = curriedFunction(tmpthis[i]); return arr; } // Otherwise run on the vector. return curriedFunction(this[0]); }; })(funcs[i]); })('quantiles percentileOfScore'.split(' ')); }(this.jStat, Math)); // Special functions // (function(jStat, Math) { // Log-gamma function jStat.gammaln = function gammaln(x) { var j = 0; var cof = [ 76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5 ]; var ser = 1.000000000190015; var xx, y, tmp; tmp = (y = xx = x) + 5.5; tmp -= (xx + 0.5) * Math.log(tmp); for (; j < 6; j++) ser += cof[j] / ++y; return Math.log(2.5066282746310005 * ser / xx) - tmp; }; // gamma of x jStat.gammafn = function gammafn(x) { var p = [-1.716185138865495, 24.76565080557592, -379.80425647094563, 629.3311553128184, 866.9662027904133, -31451.272968848367, -36144.413418691176, 66456.14382024054 ]; var q = [-30.8402300119739, 315.35062697960416, -1015.1563674902192, -3107.771671572311, 22538.118420980151, 4755.8462775278811, -134659.9598649693, -115132.2596755535]; var fact = false; var n = 0; var xden = 0; var xnum = 0; var y = x; var i, z, yi, res, sum, ysq; if (y <= 0) { res = y % 1 + 3.6e-16; if (res) { fact = (!(y & 1) ? 1 : -1) * Math.PI / Math.sin(Math.PI * res); y = 1 - y; } else { return Infinity; } } yi = y; if (y < 1) { z = y++; } else { z = (y -= n = (y | 0) - 1) - 1; } for (i = 0; i < 8; ++i) { xnum = (xnum + p[i]) * z; xden = xden * z + q[i]; } res = xnum / xden + 1; if (yi < y) { res /= yi; } else if (yi > y) { for (i = 0; i < n; ++i) { res *= y; y++; } } if (fact) { res = fact / res; } return res; }; // lower incomplete gamma function, which is usually typeset with a // lower-case greek gamma as the function symbol jStat.gammap = function gammap(a, x) { return jStat.lowRegGamma(a, x) * jStat.gammafn(a); }; // The lower regularized incomplete gamma function, usually written P(a,x) jStat.lowRegGamma = function lowRegGamma(a, x) { var aln = jStat.gammaln(a); var ap = a; var sum = 1 / a; var del = sum; var b = x + 1 - a; var c = 1 / 1.0e-30; var d = 1 / b; var h = d; var i = 1; // calculate maximum number of itterations required for a var ITMAX = -~(Math.log((a >= 1) ? a : 1 / a) * 8.5 + a * 0.4 + 17); var an, endval; if (x < 0 || a <= 0) { return NaN; } else if (x < a + 1) { for (; i <= ITMAX; i++) { sum += del *= x / ++ap; } return (sum * Math.exp(-x + a * Math.log(x) - (aln))); } for (; i <= ITMAX; i++) { an = -i * (i - a); b += 2; d = an * d + b; c = b + an / c; d = 1 / d; h *= d * c; } return (1 - h * Math.exp(-x + a * Math.log(x) - (aln))); }; // natural log factorial of n jStat.factorialln = function factorialln(n) { return n < 0 ? NaN : jStat.gammaln(n + 1); }; // factorial of n jStat.factorial = function factorial(n) { return n < 0 ? NaN : jStat.gammafn(n + 1); }; // combinations of n, m jStat.combination = function combination(n, m) { // make sure n or m don't exceed the upper limit of usable values return (n > 170 || m > 170) ? Math.exp(jStat.combinationln(n, m)) : (jStat.factorial(n) / jStat.factorial(m)) / jStat.factorial(n - m); }; jStat.combinationln = function combinationln(n, m){ return jStat.factorialln(n) - jStat.factorialln(m) - jStat.factorialln(n - m); }; // permutations of n, m jStat.permutation = function permutation(n, m) { return jStat.factorial(n) / jStat.factorial(n - m); }; // beta function jStat.betafn = function betafn(x, y) { // ensure arguments are positive if (x <= 0 || y <= 0) return undefined; // make sure x + y doesn't exceed the upper limit of usable values return (x + y > 170) ? Math.exp(jStat.betaln(x, y)) : jStat.gammafn(x) * jStat.gammafn(y) / jStat.gammafn(x + y); }; // natural logarithm of beta function jStat.betaln = function betaln(x, y) { return jStat.gammaln(x) + jStat.gammaln(y) - jStat.gammaln(x + y); }; // Evaluates the continued fraction for incomplete beta function by modified // Lentz's method. jStat.betacf = function betacf(x, a, b) { var fpmin = 1e-30; var m = 1; var qab = a + b; var qap = a + 1; var qam = a - 1; var c = 1; var d = 1 - qab * x / qap; var m2, aa, del, h; // These q's will be used in factors that occur in the coefficients if (Math.abs(d) < fpmin) d = fpmin; d = 1 / d; h = d; for (; m <= 100; m++) { m2 = 2 * m; aa = m * (b - m) * x / ((qam + m2) * (a + m2)); // One step (the even one) of the recurrence d = 1 + aa * d; if (Math.abs(d) < fpmin) d = fpmin; c = 1 + aa / c; if (Math.abs(c) < fpmin) c = fpmin; d = 1 / d; h *= d * c; aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2)); // Next step of the recurrence (the odd one) d = 1 + aa * d; if (Math.abs(d) < fpmin) d = fpmin; c = 1 + aa / c; if (Math.abs(c) < fpmin) c = fpmin; d = 1 / d; del = d * c; h *= del; if (Math.abs(del - 1.0) < 3e-7) break; } return h; }; // Returns the inverse of the lower regularized inomplete gamma function jStat.gammapinv = function gammapinv(p, a) { var j = 0; var a1 = a - 1; var EPS = 1e-8; var gln = jStat.gammaln(a); var x, err, t, u, pp, lna1, afac; if (p >= 1) return Math.max(100, a + 100 * Math.sqrt(a)); if (p <= 0) return 0; if (a > 1) { lna1 = Math.log(a1); afac = Math.exp(a1 * (lna1 - 1) - gln); pp = (p < 0.5) ? p : 1 - p; t = Math.sqrt(-2 * Math.log(pp)); x = (2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t; if (p < 0.5) x = -x; x = Math.max(1e-3, a * Math.pow(1 - 1 / (9 * a) - x / (3 * Math.sqrt(a)), 3)); } else { t = 1 - a * (0.253 + a * 0.12); if (p < t) x = Math.pow(p / t, 1 / a); else x = 1 - Math.log(1 - (p - t) / (1 - t)); } for(; j < 12; j++) { if (x <= 0) return 0; err = jStat.lowRegGamma(a, x) - p; if (a > 1) t = afac * Math.exp(-(x - a1) + a1 * (Math.log(x) - lna1)); else t = Math.exp(-x + a1 * Math.log(x) - gln); u = err / t; x -= (t = u / (1 - 0.5 * Math.min(1, u * ((a - 1) / x - 1)))); if (x <= 0) x = 0.5 * (x + t); if (Math.abs(t) < EPS * x) break; } return x; }; // Returns the error function erf(x) jStat.erf = function erf(x) { var cof = [-1.3026537197817094, 6.4196979235649026e-1, 1.9476473204185836e-2, -9.561514786808631e-3, -9.46595344482036e-4, 3.66839497852761e-4, 4.2523324806907e-5, -2.0278578112534e-5, -1.624290004647e-6, 1.303655835580e-6, 1.5626441722e-8, -8.5238095915e-8, 6.529054439e-9, 5.059343495e-9, -9.91364156e-10, -2.27365122e-10, 9.6467911e-11, 2.394038e-12, -6.886027e-12, 8.94487e-13, 3.13092e-13, -1.12708e-13, 3.81e-16, 7.106e-15, -1.523e-15, -9.4e-17, 1.21e-16, -2.8e-17]; var j = cof.length - 1; var isneg = false; var d = 0; var dd = 0; var t, ty, tmp, res; if (x < 0) { x = -x; isneg = true; } t = 2 / (2 + x); ty = 4 * t - 2; for(; j > 0; j--) { tmp = d; d = ty * d - dd + cof[j]; dd = tmp; } res = t * Math.exp(-x * x + 0.5 * (cof[0] + ty * d) - dd); return isneg ? res - 1 : 1 - res; }; // Returns the complmentary error function erfc(x) jStat.erfc = function erfc(x) { return 1 - jStat.erf(x); }; // Returns the inverse of the complementary error function jStat.erfcinv = function erfcinv(p) { var j = 0; var x, err, t, pp; if (p >= 2) return -100; if (p <= 0) return 100; pp = (p < 1) ? p : 2 - p; t = Math.sqrt(-2 * Math.log(pp / 2)); x = -0.70711 * ((2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t); for (; j < 2; j++) { err = jStat.erfc(x) - pp; x += err / (1.12837916709551257 * Math.exp(-x * x) - x * err); } return (p < 1) ? x : -x; }; // Returns the inverse of the incomplete beta function jStat.ibetainv = function ibetainv(p, a, b) { var EPS = 1e-8; var a1 = a - 1; var b1 = b - 1; var j = 0; var lna, lnb, pp, t, u, err, x, al, h, w, afac; if (p <= 0) return 0; if (p >= 1) return 1; if (a >= 1 && b >= 1) { pp = (p < 0.5) ? p : 1 - p; t = Math.sqrt(-2 * Math.log(pp)); x = (2.30753 + t * 0.27061) / (1 + t* (0.99229 + t * 0.04481)) - t; if (p < 0.5) x = -x; al = (x * x - 3) / 6; h = 2 / (1 / (2 * a - 1) + 1 / (2 * b - 1)); w = (x * Math.sqrt(al + h) / h) - (1 / (2 * b - 1) - 1 / (2 * a - 1)) * (al + 5 / 6 - 2 / (3 * h)); x = a / (a + b * Math.exp(2 * w)); } else { lna = Math.log(a / (a + b)); lnb = Math.log(b / (a + b)); t = Math.exp(a * lna) / a; u = Math.exp(b * lnb) / b; w = t + u; if (p < t / w) x = Math.pow(a * w * p, 1 / a); else x = 1 - Math.pow(b * w * (1 - p), 1 / b); } afac = -jStat.gammaln(a) - jStat.gammaln(b) + jStat.gammaln(a + b); for(; j < 10; j++) { if (x === 0 || x === 1) return x; err = jStat.ibeta(x, a, b) - p; t = Math.exp(a1 * Math.log(x) + b1 * Math.log(1 - x) + afac); u = err / t; x -= (t = u / (1 - 0.5 * Math.min(1, u * (a1 / x - b1 / (1 - x))))); if (x <= 0) x = 0.5 * (x + t); if (x >= 1) x = 0.5 * (x + t + 1); if (Math.abs(t) < EPS * x && j > 0) break; } return x; }; // Returns the incomplete beta function I_x(a,b) jStat.ibeta = function ibeta(x, a, b) { // Factors in front of the continued fraction. var bt = (x === 0 || x === 1) ? 0 : Math.exp(jStat.gammaln(a + b) - jStat.gammaln(a) - jStat.gammaln(b) + a * Math.log(x) + b * Math.log(1 - x)); if (x < 0 || x > 1) return false; if (x < (a + 1) / (a + b + 2)) // Use continued fraction directly. return bt * jStat.betacf(x, a, b) / a; // else use continued fraction after making the symmetry transformation. return 1 - bt * jStat.betacf(1 - x, b, a) / b; }; // Returns a normal deviate (mu=0, sigma=1). // If n and m are specified it returns a object of normal deviates. jStat.randn = function randn(n, m) { var u, v, x, y, q, mat; if (!m) m = n; if (n) return jStat.create(n, m, function() { return jStat.randn(); }); do { u = Math.random(); v = 1.7156 * (Math.random() - 0.5); x = u - 0.449871; y = Math.abs(v) + 0.386595; q = x * x + y * (0.19600 * y - 0.25472 * x); } while (q > 0.27597 && (q > 0.27846 || v * v > -4 * Math.log(u) * u * u)); return v / u; }; // Returns a gamma deviate by the method of Marsaglia and Tsang. jStat.randg = function randg(shape, n, m) { var oalph = shape; var a1, a2, u, v, x, mat; if (!m) m = n; if (!shape) shape = 1; if (n) { mat = jStat.zeros(n,m); mat.alter(function() { return jStat.randg(shape); }); return mat; } if (shape < 1) shape += 1; a1 = shape - 1 / 3; a2 = 1 / Math.sqrt(9 * a1); do { do { x = jStat.randn(); v = 1 + a2 * x; } while(v <= 0); v = v * v * v; u = Math.random(); } while(u > 1 - 0.331 * Math.pow(x, 4) && Math.log(u) > 0.5 * x*x + a1 * (1 - v + Math.log(v))); // alpha > 1 if (shape == oalph) return a1 * v; // alpha < 1 do { u = Math.random(); } while(u === 0); return Math.pow(u, 1 / oalph) * a1 * v; }; // making use of static methods on the instance (function(funcs) { for (var i = 0; i < funcs.length; i++) (function(passfunc) { jStat.fn[passfunc] = function() { return jStat( jStat.map(this, function(value) { return jStat[passfunc](value); })); } })(funcs[i]); })('gammaln gammafn factorial factorialln'.split(' ')); (function(funcs) { for (var i = 0; i < funcs.length; i++) (function(passfunc) { jStat.fn[passfunc] = function() { return jStat(jStat[passfunc].apply(null, arguments)); }; })(funcs[i]); })('randn'.split(' ')); }(this.jStat, Math)); (function(jStat, Math) { // generate all distribution instance methods (function(list) { for (var i = 0; i < list.length; i++) (function(func) { // distribution instance method jStat[func] = function(a, b, c) { if (!(this instanceof arguments.callee)) return new arguments.callee(a, b, c); this._a = a; this._b = b; this._c = c; return this; }; // distribution method to be used on a jStat instance jStat.fn[func] = function(a, b, c) { var newthis = jStat[func](a, b, c); newthis.data = this; return newthis; }; // sample instance method jStat[func].prototype.sample = function(arr) { var a = this._a; var b = this._b; var c = this._c; if (arr) return jStat.alter(arr, function() { return jStat[func].sample(a, b, c); }); else return jStat[func].sample(a, b, c); }; // generate the pdf, cdf and inv instance methods (function(vals) { for (var i = 0; i < vals.length; i++) (function(fnfunc) { jStat[func].prototype[fnfunc] = function(x) { var a = this._a; var b = this._b; var c = this._c; if (!x && x !== 0) x = this.data; if (typeof x !== 'number') { return jStat.fn.map.call(x, function(x) { return jStat[func][fnfunc](x, a, b, c); }); } return jStat[func][fnfunc](x, a, b, c); }; })(vals[i]); })('pdf cdf inv'.split(' ')); // generate the mean, median, mode and variance instance methods (function(vals) { for (var i = 0; i < vals.length; i++) (function(fnfunc) { jStat[func].prototype[fnfunc] = function() { return jStat[func][fnfunc](this._a, this._b, this._c); }; })(vals[i]); })('mean median mode variance'.split(' ')); })(list[i]); })(( 'beta centralF cauchy chisquare exponential gamma invgamma kumaraswamy ' + 'lognormal normal pareto studentt weibull uniform binomial negbin hypgeom ' + 'poisson triangular' ).split(' ')); // extend beta function with static methods jStat.extend(jStat.beta, { pdf: function pdf(x, alpha, beta) { // PDF is zero outside the support if (x > 1 || x < 0) return 0; // PDF is one for the uniform case if (alpha == 1 && beta == 1) return 1; if (alpha < 512 || beta < 512) { return (Math.pow(x, alpha - 1) * Math.pow(1 - x, beta - 1)) / jStat.betafn(alpha, beta); } else { return Math.exp((alpha - 1) * Math.log(x) + (beta - 1) * Math.log(1 - x) - jStat.betaln(alpha, beta)); } }, cdf: function cdf(x, alpha, beta) { return (x > 1 || x < 0) ? (x > 1) * 1 : jStat.ibeta(x, alpha, beta); }, inv: function inv(x, alpha, beta) { return jStat.ibetainv(x, alpha, beta); }, mean: function mean(alpha, beta) { return alpha / (alpha + beta); }, median: function median(alpha, beta) { throw new Error('median not yet implemented'); }, mode: function mode(alpha, beta) { return (alpha * beta) / (Math.pow(alpha + beta, 2) * (alpha + beta + 1)); }, // return a random sample sample: function sample(alpha, beta) { var u = jStat.randg(alpha); return u / (u + jStat.randg(beta)); }, variance: function variance(alpha, beta) { return (alpha * beta) / (Math.pow(alpha + beta, 2) * (alpha + beta + 1)); } }); // extend F function with static methods jStat.extend(jStat.centralF, { // This implementation of the pdf function avoids float overflow // See the way that R calculates this value: // https://svn.r-project.org/R/trunk/src/nmath/df.c pdf: function pdf(x, df1, df2) { var p, q, f; if (x < 0) return undefined; if (df1 <= 2) { if (df1 === 1 && df2 === 1) { return Infinity; } if (df1 === 2 && df2 === 1) { return 1; } return Math.sqrt((Math.pow(df1 * x, df1) * Math.pow(df2, df2)) / (Math.pow(df1 * x + df2, df1 + df2))) / (x * jStat.betafn(df1/2, df2/2)); } p = (df1 * x) / (df2 + x * df1); q = df2 / (df2 + x * df1); f = df1 * q / 2.0; return f * jStat.binomial.pdf((df1 - 2) / 2, (df1 + df2 - 2) / 2, p); }, cdf: function cdf(x, df1, df2) { return jStat.ibeta((df1 * x) / (df1 * x + df2), df1 / 2, df2 / 2); }, inv: function inv(x, df1, df2) { return df2 / (df1 * (1 / jStat.ibetainv(x, df1 / 2, df2 / 2) - 1)); }, mean: function mean(df1, df2) { return (df2 > 2) ? df2 / (df2 - 2) : undefined; }, mode: function mode(df1, df2) { return (df1 > 2) ? (df2 * (df1 - 2)) / (df1 * (df2 + 2)) : undefined; }, // return a random sample sample: function sample(df1, df2) { var x1 = jStat.randg(df1 / 2) * 2; var x2 = jStat.randg(df2 / 2) * 2; return (x1 / df1) / (x2 / df2); }, variance: function variance(df1, df2) { if (df2 <= 4) return undefined; return 2 * df2 * df2 * (df1 + df2 - 2) / (df1 * (df2 - 2) * (df2 - 2) * (df2 - 4)); } }); // extend cauchy function with static methods jStat.extend(jStat.cauchy, { pdf: function pdf(x, local, scale) { return (scale / (Math.pow(x - local, 2) + Math.pow(scale, 2))) / Math.PI; }, cdf: function cdf(x, local, scale) { return Math.atan((x - local) / scale) / Math.PI + 0.5; }, inv: function(p, local, scale) { return local + scale * Math.tan(Math.PI * (p - 0.5)); }, median: function median(local, scale) { return local; }, mode: function mode(local, scale) { return local; }, sample: function sample(local, scale) { return jStat.randn() * Math.sqrt(1 / (2 * jStat.randg(0.5))) * scale + local; } }); // extend chisquare function with static methods jStat.extend(jStat.chisquare, { pdf: function pdf(x, dof) { return x === 0 ? 0 : Math.exp((dof / 2 - 1) * Math.log(x) - x / 2 - (dof / 2) * Math.log(2) - jStat.gammaln(dof / 2)); }, cdf: function cdf(x, dof) { return jStat.lowRegGamma(dof / 2, x / 2); }, inv: function(p, dof) { return 2 * jStat.gammapinv(p, 0.5 * dof); }, mean : function(dof) { return dof; }, // TODO: this is an approximation (is there a better way?) median: function median(dof) { return dof * Math.pow(1 - (2 / (9 * dof)), 3); }, mode: function mode(dof) { return (dof - 2 > 0) ? dof - 2 : 0; }, sample: function sample(dof) { return jStat.randg(dof / 2) * 2; }, variance: function variance(dof) { return 2 * dof; } }); // extend exponential function with static methods jStat.extend(jStat.exponential, { pdf: function pdf(x, rate) { return x < 0 ? 0 : rate * Math.exp(-rate * x); }, cdf: function cdf(x, rate) { return x < 0 ? 0 : 1 - Math.exp(-rate * x); }, inv: function(p, rate) { return -Math.log(1 - p) / rate; }, mean : function(rate) { return 1 / rate; }, median: function (rate) { return (1 / rate) * Math.log(2); }, mode: function mode(rate) { return 0; }, sample: function sample(rate) { return -1 / rate * Math.log(Math.random()); }, variance : function(rate) { return Math.pow(rate, -2); } }); // extend gamma function with static methods jStat.extend(jStat.gamma, { pdf: function pdf(x, shape, scale) { return Math.exp((shape - 1) * Math.log(x) - x / scale - jStat.gammaln(shape) - shape * Math.log(scale)); }, cdf: function cdf(x, shape, scale) { return jStat.lowRegGamma(shape, x / scale); }, inv: function(p, shape, scale) { return jStat.gammapinv(p, shape) * scale; }, mean : function(shape, scale) { return shape * scale; }, mode: function mode(shape, scale) { if(shape > 1) return (shape - 1) * scale; return undefined; }, sample: function sample(shape, scale) { return jStat.randg(shape) * scale; }, variance: function variance(shape, scale) { return shape * scale * scale; } }); // extend inverse gamma function with static methods jStat.extend(jStat.invgamma, { pdf: function pdf(x, shape, scale) { return Math.exp(-(shape + 1) * Math.log(x) - scale / x - jStat.gammaln(shape) + shape * Math.log(scale)); }, cdf: function cdf(x, shape, scale) { return 1 - jStat.lowRegGamma(shape, scale / x); }, inv: function(p, shape, scale) { return scale / jStat.gammapinv(1 - p, shape); }, mean : function(shape, scale) { return (shape > 1) ? scale / (shape - 1) : undefined; }, mode: function mode(shape, scale) { return scale / (shape + 1); }, sample: function sample(shape, scale) { return scale / jStat.randg(shape); }, variance: function variance(shape, scale) { if (shape <= 2) return undefined; return scale * scale / ((shape - 1) * (shape - 1) * (shape - 2)); } }); // extend kumaraswamy function with static methods jStat.extend(jStat.kumaraswamy, { pdf: function pdf(x, alpha, beta) { return Math.exp(Math.log(alpha) + Math.log(beta) + (alpha - 1) * Math.log(x) + (beta - 1) * Math.log(1 - Math.pow(x, alpha))); }, cdf: function cdf(x, alpha, beta) { return (1 - Math.pow(1 - Math.pow(x, alpha), beta)); }, mean : function(alpha, beta) { return (beta * jStat.gammafn(1 + 1 / alpha) * jStat.gammafn(beta)) / (jStat.gammafn(1 + 1 / alpha + beta)); }, median: function median(alpha, beta) { return Math.pow(1 - Math.pow(2, -1 / beta), 1 / alpha); }, mode: function mode(alpha, beta) { if (!(alpha >= 1 && beta >= 1 && (alpha !== 1 && beta !== 1))) return undefined; return Math.pow((alpha - 1) / (alpha * beta - 1), 1 / alpha); }, variance: function variance(alpha, beta) { throw new Error('variance not yet implemented'); // TODO: complete this } }); // extend lognormal function with static methods jStat.extend(jStat.lognormal, { pdf: function pdf(x, mu, sigma) { return Math.exp(-Math.log(x) - 0.5 * Math.log(2 * Math.PI) - Math.log(sigma) - Math.pow(Math.log(x) - mu, 2) / (2 * sigma * sigma)); }, cdf: function cdf(x, mu, sigma) { return 0.5 + (0.5 * jStat.erf((Math.log(x) - mu) / Math.sqrt(2 * sigma * sigma))); }, inv: function(p, mu, sigma) { return Math.exp(-1.41421356237309505 * sigma * jStat.erfcinv(2 * p) + mu); }, mean: function mean(mu, sigma) { return Math.exp(mu + sigma * sigma / 2); }, median: function median(mu, sigma) { return Math.exp(mu); }, mode: function mode(mu, sigma) { return Math.exp(mu - sigma * sigma); }, sample: function sample(mu, sigma) { return Math.exp(jStat.randn() * sigma + mu); }, variance: function variance(mu, sigma) { return (Math.exp(sigma * sigma) - 1) * Math.exp(2 * mu + sigma * sigma); } }); // extend normal function with static methods jStat.extend(jStat.normal, { pdf: function pdf(x, mean, std) { return Math.exp(-0.5 * Math.log(2 * Math.PI) - Math.log(std) - Math.pow(x - mean, 2) / (2 * std * std)); }, cdf: function cdf(x, mean, std) { return 0.5 * (1 + jStat.erf((x - mean) / Math.sqrt(2 * std * std))); }, inv: function(p, mean, std) { return -1.41421356237309505 * std * jStat.erfcinv(2 * p) + mean; }, mean : function(mean, std) { return mean; }, median: function median(mean, std) { return mean; }, mode: function (mean, std) { return mean; }, sample: function sample(mean, std) { return jStat.randn() * std + mean; }, variance : function(mean, std) { return std * std; } }); // extend pareto function with static methods jStat.extend(jStat.pareto, { pdf: function pdf(x, scale, shape) { if (x < scale) return undefined; return (shape * Math.pow(scale, shape)) / Math.pow(x, shape + 1); }, cdf: function cdf(x, scale, shape) { return 1 - Math.pow(scale / x, shape); }, mean: function mean(scale, shape) { if (shape <= 1) return undefined; return (shape * Math.pow(scale, shape)) / (shape - 1); }, median: function median(scale, shape) { return scale * (shape * Math.SQRT2); }, mode: function mode(scale, shape) { return scale; }, variance : function(scale, shape) { if (shape <= 2) return undefined; return (scale*scale * shape) / (Math.pow(shape - 1, 2) * (shape - 2)); } }); // extend studentt function with static methods jStat.extend(jStat.studentt, { pdf: function pdf(x, dof) { return (jStat.gammafn((dof + 1) / 2) / (Math.sqrt(dof * Math.PI) * jStat.gammafn(dof / 2))) * Math.pow(1 + ((x * x) / dof), -((dof + 1) / 2)); }, cdf: function cdf(x, dof) { var dof2 = dof / 2; return jStat.ibeta((x + Math.sqrt(x * x + dof)) / (2 * Math.sqrt(x * x + dof)), dof2, dof2); }, inv: function(p, dof) { var x = jStat.ibetainv(2 * Math.min(p, 1 - p), 0.5 * dof, 0.5); x = Math.sqrt(dof * (1 - x) / x); return (p > 0.5) ? x : -x; }, mean: function mean(dof) { return (dof > 1) ? 0 : undefined; }, median: function median(dof) { return 0; }, mode: function mode(dof) { return 0; }, sample: function sample(dof) { return jStat.randn() * Math.sqrt(dof / (2 * jStat.randg(dof / 2))); }, variance: function variance(dof) { return (dof > 2) ? dof / (dof - 2) : (dof > 1) ? Infinity : undefined; } }); // extend weibull function with static methods jStat.extend(jStat.weibull, { pdf: function pdf(x, scale, shape) { if (x < 0) return 0; return (shape / scale) * Math.pow((x / scale), (shape - 1)) * Math.exp(-(Math.pow((x / scale), shape))); }, cdf: function cdf(x, scale, shape) { return x < 0 ? 0 : 1 - Math.exp(-Math.pow((x / scale), shape)); }, inv: function(p, scale, shape) { return scale * Math.pow(-Math.log(1 - p), 1 / shape); }, mean : function(scale, shape) { return scale * jStat.gammafn(1 + 1 / shape); }, median: function median(scale, shape) { return scale * Math.pow(Math.log(2), 1 / shape); }, mode: function mode(scale, shape) { if (shape <= 1) return undefined; return scale * Math.pow((shape - 1) / shape, 1 / shape); }, sample: function sample(scale, shape) { return scale * Math.pow(-Math.log(Math.random()), 1 / shape); }, variance: function variance(scale, shape) { return scale * scale * jStat.gammafn(1 + 2 / shape) - Math.pow(this.mean(scale, shape), 2); } }); // extend uniform function with static methods jStat.extend(jStat.uniform, { pdf: function pdf(x, a, b) { return (x < a || x > b) ? 0 : 1 / (b - a); }, cdf: function cdf(x, a, b) { if (x < a) return 0; else if (x < b) return (x - a) / (b - a); return 1; }, inv: function(p, a, b) { return a + (p * (b - a)); }, mean: function mean(a, b) { return 0.5 * (a + b); }, median: function median(a, b) { return jStat.mean(a, b); }, mode: function mode(a, b) { throw new Error('mode is not yet implemented'); }, sample: function sample(a, b) { return (a / 2 + b / 2) + (b / 2 - a / 2) * (2 * Math.random() - 1); }, variance: function variance(a, b) { return Math.pow(b - a, 2) / 12; } }); // extend uniform function with static methods jStat.extend(jStat.binomial, { pdf: function pdf(k, n, p) { return (p === 0 || p === 1) ? ((n * p) === k ? 1 : 0) : jStat.combination(n, k) * Math.pow(p, k) * Math.pow(1 - p, n - k); }, cdf: function cdf(x, n, p) { var binomarr = [], k = 0; if (x < 0) { return 0; } if (x < n) { for (; k <= x; k++) { binomarr[ k ] = jStat.binomial.pdf(k, n, p); } return jStat.sum(binomarr); } return 1; } }); // extend uniform function with static methods jStat.extend(jStat.negbin, { pdf: function pdf(k, r, p) { return k !== k | 0 ? false : k < 0 ? 0 : jStat.combination(k + r - 1, r - 1) * Math.pow(1 - p, k) * Math.pow(p, r); }, cdf: function cdf(x, r, p) { var sum = 0, k = 0; if (x < 0) return 0; for (; k <= x; k++) { sum += jStat.negbin.pdf(k, r, p); } return sum; } }); // extend uniform function with static methods jStat.extend(jStat.hypgeom, { pdf: function pdf(k, N, m, n) { // Hypergeometric PDF. // A simplification of the CDF algorithm below. // k = number of successes drawn // N = population size // m = number of successes in population // n = number of items drawn from population if(k !== k | 0) { return false; } else if(k < 0 || k < m - (N - n)) { // It's impossible to have this few successes drawn. return 0; } else if(k > n || k > m) { // It's impossible to have this many successes drawn. return 0; } else if (m * 2 > N) { // More than half the population is successes. if(n * 2 > N) { // More than half the population is sampled. return jStat.hypgeom.pdf(N - m - n + k, N, N - m, N - n) } else { // Half or less of the population is sampled. return jStat.hypgeom.pdf(n - k, N, N - m, n); } } else if(n * 2 > N) { // Half or less is successes. return jStat.hypgeom.pdf(m - k, N, m, N - n); } else if(m < n) { // We want to have the number of things sampled to be less than the // successes available. So swap the definitions of successful and sampled. return jStat.hypgeom.pdf(k, N, n, m); } else { // If we get here, half or less of the population was sampled, half or // less of it was successes, and we had fewer sampled things than // successes. Now we can do this complicated iterative algorithm in an // efficient way. // The basic premise of the algorithm is that we partially normalize our // intermediate product to keep it in a numerically good region, and then // finish the normalization at the end. // This variable holds the scaled probability of the current number of // successes. var scaledPDF = 1; // This keeps track of how much we have normalized. var samplesDone = 0; for(var i = 0; i < k; i++) { // For every possible number of successes up to that observed... while(scaledPDF > 1 && samplesDone < n) { // Intermediate result is growing too big. Apply some of the // normalization to shrink everything. scaledPDF *= 1 - (m / (N - samplesDone)); // Say we've normalized by this sample already. samplesDone++; } // Work out the partially-normalized hypergeometric PDF for the next // number of successes scaledPDF *= (n - i) * (m - i) / ((i + 1) * (N - m - n + i + 1)); } for(; samplesDone < n; samplesDone++) { // Apply all the rest of the normalization scaledPDF *= 1 - (m / (N - samplesDone)); } // Bound answer sanely before returning. return Math.min(1, Math.max(0, scaledPDF)); } }, cdf: function cdf(x, N, m, n) { // Hypergeometric CDF. // This algorithm is due to Prof. Thomas S. Ferguson, , // and comes from his hypergeometric test calculator at // . // x = number of successes drawn // N = population size // m = number of successes in population // n = number of items drawn from population if(x < 0 || x < m - (N - n)) { // It's impossible to have this few successes drawn or fewer. return 0; } else if(x >= n || x >= m) { // We will always have this many successes or fewer. return 1; } else if (m * 2 > N) { // More than half the population is successes. if(n * 2 > N) { // More than half the population is sampled. return jStat.hypgeom.cdf(N - m - n + x, N, N - m, N - n) } else { // Half or less of the population is sampled. return 1 - jStat.hypgeom.cdf(n - x - 1, N, N - m, n); } } else if(n * 2 > N) { // Half or less is successes. return 1 - jStat.hypgeom.cdf(m - x - 1, N, m, N - n); } else if(m < n) { // We want to have the number of things sampled to be less than the // successes available. So swap the definitions of successful and sampled. return jStat.hypgeom.cdf(x, N, n, m); } else { // If we get here, half or less of the population was sampled, half or // less of it was successes, and we had fewer sampled things than // successes. Now we can do this complicated iterative algorithm in an // efficient way. // The basic premise of the algorithm is that we partially normalize our // intermediate sum to keep it in a numerically good region, and then // finish the normalization at the end. // Holds the intermediate, scaled total CDF. var scaledCDF = 1; // This variable holds the scaled probability of the current number of // successes. var scaledPDF = 1; // This keeps track of how much we have normalized. var samplesDone = 0; for(var i = 0; i < x; i++) { // For every possible number of successes up to that observed... while(scaledCDF > 1 && samplesDone < n) { // Intermediate result is growing too big. Apply some of the // normalization to shrink everything. var factor = 1 - (m / (N - samplesDone)); scaledPDF *= factor; scaledCDF *= factor; // Say we've normalized by this sample already. samplesDone++; } // Work out the partially-normalized hypergeometric PDF for the next // number of successes scaledPDF *= (n - i) * (m - i) / ((i + 1) * (N - m - n + i + 1)); // Add to the CDF answer. scaledCDF += scaledPDF; } for(; samplesDone < n; samplesDone++) { // Apply all the rest of the normalization scaledCDF *= 1 - (m / (N - samplesDone)); } // Bound answer sanely before returning. return Math.min(1, Math.max(0, scaledCDF)); } } }); // extend uniform function with static methods jStat.extend(jStat.poisson, { pdf: function pdf(k, l) { return Math.pow(l, k) * Math.exp(-l) / jStat.factorial(k); }, cdf: function cdf(x, l) { var sumarr = [], k = 0; if (x < 0) return 0; for (; k <= x; k++) { sumarr.push(jStat.poisson.pdf(k, l)); } return jStat.sum(sumarr); }, mean : function(l) { return l; }, variance : function(l) { return l; }, sample: function sample(l) { var p = 1, k = 0, L = Math.exp(-l); do { k++; p *= Math.random(); } while (p > L); return k - 1; } }); // extend triangular function with static methods jStat.extend(jStat.triangular, { pdf: function pdf(x, a, b, c) { return (b <= a || c < a || c > b) ? undefined : (x < a || x > b) ? 0 : (x <= c) ? (2 * (x - a)) / ((b - a) * (c - a)) : (2 * (b - x)) / ((b - a) * (b - c)); }, cdf: function cdf(x, a, b, c) { if (b <= a || c < a || c > b) return undefined; if (x < a) { return 0; } else { if (x <= c) return Math.pow(x - a, 2) / ((b - a) * (c - a)); return 1 - Math.pow(b - x, 2) / ((b - a) * (b - c)); } // never reach this return 1; }, mean: function mean(a, b, c) { return (a + b + c) / 3; }, median: function median(a, b, c) { if (c <= (a + b) / 2) { return b - Math.sqrt((b - a) * (b - c)) / Math.sqrt(2); } else if (c > (a + b) / 2) { return a + Math.sqrt((b - a) * (c - a)) / Math.sqrt(2); } }, mode: function mode(a, b, c) { return c; }, sample: function sample(a, b, c) { var u = Math.random(); if (u < ((c - a) / (b - a))) return a + Math.sqrt(u * (b - a) * (c - a)) return b - Math.sqrt((1 - u) * (b - a) * (b - c)); }, variance: function variance(a, b, c) { return (a * a + b * b + c * c - a * b - a * c - b * c) / 18; } }); }(this.jStat, Math)); /* Provides functions for the solution of linear system of equations, integration, extrapolation, * interpolation, eigenvalue problems, differential equations and PCA analysis. */ (function(jStat, Math) { var push = Array.prototype.push; var isArray = jStat.utils.isArray; jStat.extend({ // add a vector/matrix to a vector/matrix or scalar add: function add(arr, arg) { // check if arg is a vector or scalar if (isArray(arg)) { if (!isArray(arg[0])) arg = [ arg ]; return jStat.map(arr, function(value, row, col) { return value + arg[row][col]; }); } return jStat.map(arr, function(value) { return value + arg; }); }, // subtract a vector or scalar from the vector subtract: function subtract(arr, arg) { // check if arg is a vector or scalar if (isArray(arg)) { if (!isArray(arg[0])) arg = [ arg ]; return jStat.map(arr, function(value, row, col) { return value - arg[row][col] || 0; }); } return jStat.map(arr, function(value) { return value - arg; }); }, // matrix division divide: function divide(arr, arg) { if (isArray(arg)) { if (!isArray(arg[0])) arg = [ arg ]; return jStat.multiply(arr, jStat.inv(arg)); } return jStat.map(arr, function(value) { return value / arg; }); }, // matrix multiplication multiply: function multiply(arr, arg) { var row, col, nrescols, sum, nrow = arr.length, ncol = arr[0].length, res = jStat.zeros(nrow, nrescols = (isArray(arg)) ? arg[0].length : ncol), rescols = 0; if (isArray(arg)) { for (; rescols < nrescols; rescols++) { for (row = 0; row < nrow; row++) { sum = 0; for (col = 0; col < ncol; col++) sum += arr[row][col] * arg[col][rescols]; res[row][rescols] = sum; } } return (nrow === 1 && rescols === 1) ? res[0][0] : res; } return jStat.map(arr, function(value) { return value * arg; }); }, // Returns the dot product of two matricies dot: function dot(arr, arg) { if (!isArray(arr[0])) arr = [ arr ]; if (!isArray(arg[0])) arg = [ arg ]; // convert column to row vector var left = (arr[0].length === 1 && arr.length !== 1) ? jStat.transpose(arr) : arr, right = (arg[0].length === 1 && arg.length !== 1) ? jStat.transpose(arg) : arg, res = [], row = 0, nrow = left.length, ncol = left[0].length, sum, col; for (; row < nrow; row++) { res[row] = []; sum = 0; for (col = 0; col < ncol; col++) sum += left[row][col] * right[row][col]; res[row] = sum; } return (res.length === 1) ? res[0] : res; }, // raise every element by a scalar pow: function pow(arr, arg) { return jStat.map(arr, function(value) { return Math.pow(value, arg); }); }, // exponentiate every element exp: function exp(arr) { return jStat.map(arr, function(value) { return Math.exp(value); }); }, // generate the natural log of every element log: function exp(arr) { return jStat.map(arr, function(value) { return Math.log(value); }); }, // generate the absolute values of the vector abs: function abs(arr) { return jStat.map(arr, function(value) { return Math.abs(value); }); }, // computes the p-norm of the vector // In the case that a matrix is passed, uses the first row as the vector norm: function norm(arr, p) { var nnorm = 0, i = 0; // check the p-value of the norm, and set for most common case if (isNaN(p)) p = 2; // check if multi-dimensional array, and make vector correction if (isArray(arr[0])) arr = arr[0]; // vector norm for (; i < arr.length; i++) { nnorm += Math.pow(Math.abs(arr[i]), p); } return Math.pow(nnorm, 1 / p); }, // computes the angle between two vectors in rads // In case a matrix is passed, this uses the first row as the vector angle: function angle(arr, arg) { return Math.acos(jStat.dot(arr, arg) / (jStat.norm(arr) * jStat.norm(arg))); }, // augment one matrix by another // Note: this function returns a matrix, not a jStat object aug: function aug(a, b) { var newarr = a.slice(), i = 0; for (; i < newarr.length; i++) { push.apply(newarr[i], b[i]); } return newarr; }, // The inv() function calculates the inverse of a matrix // Create the inverse by augmenting the matrix by the identity matrix of the // appropriate size, and then use G-J elimination on the augmented matrix. inv: function inv(a) { var rows = a.length; var cols = a[0].length; var b = jStat.identity(rows, cols); var c = jStat.gauss_jordan(a, b); var result = []; var i = 0; var j; //We need to copy the inverse portion to a new matrix to rid G-J artifacts for (; i < rows; i++) { result[i] = []; for (j = cols; j < c[0].length; j++) result[i][j - cols] = c[i][j]; } return result; }, // calculate the determinant of a matrix det: function det(a) { var alen = a.length, alend = alen * 2, vals = new Array(alend), rowshift = alen - 1, colshift = alend - 1, mrow = rowshift - alen + 1, mcol = colshift, i = 0, result = 0, j; // check for special 2x2 case if (alen === 2) { return a[0][0] * a[1][1] - a[0][1] * a[1][0]; } for (; i < alend; i++) { vals[i] = 1; } for (i = 0; i < alen; i++) { for (j = 0; j < alen; j++) { vals[(mrow < 0) ? mrow + alen : mrow ] *= a[i][j]; vals[(mcol < alen) ? mcol + alen : mcol ] *= a[i][j]; mrow++; mcol--; } mrow = --rowshift - alen + 1; mcol = --colshift; } for (i = 0; i < alen; i++) { result += vals[i]; } for (; i < alend; i++) { result -= vals[i]; } return result; }, gauss_elimination: function gauss_elimination(a, b) { var i = 0, j = 0, n = a.length, m = a[0].length, factor = 1, sum = 0, x = [], maug, pivot, temp, k; a = jStat.aug(a, b); maug = a[0].length; for(i = 0; i < n; i++) { pivot = a[i][i]; j = i; for (k = i + 1; k < m; k++) { if (pivot < Math.abs(a[k][i])) { pivot = a[k][i]; j = k; } } if (j != i) { for(k = 0; k < maug; k++) { temp = a[i][k]; a[i][k] = a[j][k]; a[j][k] = temp; } } for (j = i + 1; j < n; j++) { factor = a[j][i] / a[i][i]; for(k = i; k < maug; k++) { a[j][k] = a[j][k] - factor * a[i][k]; } } } for (i = n - 1; i >= 0; i--) { sum = 0; for (j = i + 1; j<= n - 1; j++) { sum = sum + x[j] * a[i][j]; } x[i] =(a[i][maug - 1] - sum) / a[i][i]; } return x; }, gauss_jordan: function gauss_jordan(a, b) { var m = jStat.aug(a, b), h = m.length, w = m[0].length; // find max pivot for (var y = 0; y < h; y++) { var maxrow = y; for (var y2 = y+1; y2 < h; y2++) { if (Math.abs(m[y2][y]) > Math.abs(m[maxrow][y])) maxrow = y2; } var tmp = m[y]; m[y] = m[maxrow]; m[maxrow] = tmp for (var y2 = y+1; y2 < h; y2++) { c = m[y2][y] / m[y][y]; for (var x = y; x < w; x++) { m[y2][x] -= m[y][x] * c; } } } // backsubstitute for (var y = h-1; y >= 0; y--) { c = m[y][y]; for (var y2 = 0; y2 < y; y2++) { for (var x = w-1; x > y-1; x--) { m[y2][x] -= m[y][x] * m[y2][y] / c; } } m[y][y] /= c; for (var x = h; x < w; x++) { m[y][x] /= c; } } return m; }, lu: function lu(a, b) { throw new Error('lu not yet implemented'); }, cholesky: function cholesky(a, b) { throw new Error('cholesky not yet implemented'); }, gauss_jacobi: function gauss_jacobi(a, b, x, r) { var i = 0; var j = 0; var n = a.length; var l = []; var u = []; var d = []; var xv, c, h, xk; for (; i < n; i++) { l[i] = []; u[i] = []; d[i] = []; for (j = 0; j < n; j++) { if (i > j) { l[i][j] = a[i][j]; u[i][j] = d[i][j] = 0; } else if (i < j) { u[i][j] = a[i][j]; l[i][j] = d[i][j] = 0; } else { d[i][j] = a[i][j]; l[i][j] = u[i][j] = 0; } } } h = jStat.multiply(jStat.multiply(jStat.inv(d), jStat.add(l, u)), -1); c = jStat.multiply(jStat.inv(d), b); xv = x; xk = jStat.add(jStat.multiply(h, x), c); i = 2; while (Math.abs(jStat.norm(jStat.subtract(xk,xv))) > r) { xv = xk; xk = jStat.add(jStat.multiply(h, xv), c); i++; } return xk; }, gauss_seidel: function gauss_seidel(a, b, x, r) { var i = 0; var n = a.length; var l = []; var u = []; var d = []; var j, xv, c, h, xk; for (; i < n; i++) { l[i] = []; u[i] = []; d[i] = []; for (j = 0; j < n; j++) { if (i > j) { l[i][j] = a[i][j]; u[i][j] = d[i][j] = 0; } else if (i < j) { u[i][j] = a[i][j]; l[i][j] = d[i][j] = 0; } else { d[i][j] = a[i][j]; l[i][j] = u[i][j] = 0; } } } h = jStat.multiply(jStat.multiply(jStat.inv(jStat.add(d, l)), u), -1); c = jStat.multiply(jStat.inv(jStat.add(d, l)), b); xv = x; xk = jStat.add(jStat.multiply(h, x), c); i = 2; while (Math.abs(jStat.norm(jStat.subtract(xk, xv))) > r) { xv = xk; xk = jStat.add(jStat.multiply(h, xv), c); i = i + 1; } return xk; }, SOR: function SOR(a, b, x, r, w) { var i = 0; var n = a.length; var l = []; var u = []; var d = []; var j, xv, c, h, xk; for (; i < n; i++) { l[i] = []; u[i] = []; d[i] = []; for (j = 0; j < n; j++) { if (i > j) { l[i][j] = a[i][j]; u[i][j] = d[i][j] = 0; } else if (i < j) { u[i][j] = a[i][j]; l[i][j] = d[i][j] = 0; } else { d[i][j] = a[i][j]; l[i][j] = u[i][j] = 0; } } } h = jStat.multiply(jStat.inv(jStat.add(d, jStat.multiply(l, w))), jStat.subtract(jStat.multiply(d, 1 - w), jStat.multiply(u, w))); c = jStat.multiply(jStat.multiply(jStat.inv(jStat.add(d, jStat.multiply(l, w))), b), w); xv = x; xk = jStat.add(jStat.multiply(h, x), c); i = 2; while (Math.abs(jStat.norm(jStat.subtract(xk, xv))) > r) { xv = xk; xk = jStat.add(jStat.multiply(h, xv), c); i++; } return xk; }, householder: function householder(a) { var m = a.length; var n = a[0].length; var i = 0; var w = []; var p = []; var alpha, r, k, j, factor; for (; i < m - 1; i++) { alpha = 0; for (j = i + 1; j < n; j++) alpha += (a[j][i] * a[j][i]); factor = (a[i + 1][i] > 0) ? -1 : 1; alpha = factor * Math.sqrt(alpha); r = Math.sqrt((((alpha * alpha) - a[i + 1][i] * alpha) / 2)); w = jStat.zeros(m, 1); w[i + 1][0] = (a[i + 1][i] - alpha) / (2 * r); for (k = i + 2; k < m; k++) w[k][0] = a[k][i] / (2 * r); p = jStat.subtract(jStat.identity(m, n), jStat.multiply(jStat.multiply(w, jStat.transpose(w)), 2)); a = jStat.multiply(p, jStat.multiply(a, p)); } return a; }, // TODO: not working properly. QR: function QR(a, b) { var m = a.length; var n = a[0].length; var i = 0; var w = []; var p = []; var x = []; var j, alpha, r, k, factor, sum; for (; i < m - 1; i++) { alpha = 0; for (j = i + 1; j < n; j++) alpha += (a[j][i] * a[j][i]); factor = (a[i + 1][i] > 0) ? -1 : 1; alpha = factor * Math.sqrt(alpha); r = Math.sqrt((((alpha * alpha) - a[i + 1][i] * alpha) / 2)); w = jStat.zeros(m, 1); w[i + 1][0] = (a[i + 1][i] - alpha) / (2 * r); for (k = i + 2; k < m; k++) w[k][0] = a[k][i] / (2 * r); p = jStat.subtract(jStat.identity(m, n), jStat.multiply(jStat.multiply(w, jStat.transpose(w)), 2)); a = jStat.multiply(p, a); b = jStat.multiply(p, b); } for (i = m - 1; i >= 0; i--) { sum = 0; for (j = i + 1; j <= n - 1; j++) sum = x[j] * a[i][j]; x[i] = b[i][0] / a[i][i]; } return x; }, jacobi: function jacobi(a) { var condition = 1; var count = 0; var n = a.length; var e = jStat.identity(n, n); var ev = []; var b, i, j, p, q, maxim, theta, s; // condition === 1 only if tolerance is not reached while (condition === 1) { count++; maxim = a[0][1]; p = 0; q = 1; for (i = 0; i < n; i++) { for (j = 0; j < n; j++) { if (i != j) { if (maxim < Math.abs(a[i][j])) { maxim = Math.abs(a[i][j]); p = i; q = j; } } } } if (a[p][p] === a[q][q]) theta = (a[p][q] > 0) ? Math.PI / 4 : -Math.PI / 4; else theta = Math.atan(2 * a[p][q] / (a[p][p] - a[q][q])) / 2; s = jStat.identity(n, n); s[p][p] = Math.cos(theta); s[p][q] = -Math.sin(theta); s[q][p] = Math.sin(theta); s[q][q] = Math.cos(theta); // eigen vector matrix e = jStat.multiply(e, s); b = jStat.multiply(jStat.multiply(jStat.inv(s), a), s); a = b; condition = 0; for (i = 1; i < n; i++) { for (j = 1; j < n; j++) { if (i != j && Math.abs(a[i][j]) > 0.001) { condition = 1; } } } } for (i = 0; i < n; i++) ev.push(a[i][i]); //returns both the eigenvalue and eigenmatrix return [e, ev]; }, rungekutta: function rungekutta(f, h, p, t_j, u_j, order) { var k1, k2, u_j1, k3, k4; if (order === 2) { while (t_j <= p) { k1 = h * f(t_j, u_j); k2 = h * f(t_j + h, u_j + k1); u_j1 = u_j + (k1 + k2) / 2; u_j = u_j1; t_j = t_j + h; } } if (order === 4) { while (t_j <= p) { k1 = h * f(t_j, u_j); k2 = h * f(t_j + h / 2, u_j + k1 / 2); k3 = h * f(t_j + h / 2, u_j + k2 / 2); k4 = h * f(t_j +h, u_j + k3); u_j1 = u_j + (k1 + 2 * k2 + 2 * k3 + k4) / 6; u_j = u_j1; t_j = t_j + h; } } return u_j; }, romberg: function romberg(f, a, b, order) { var i = 0; var h = (b - a) / 2; var x = []; var h1 = []; var g = []; var m, a1, j, k, I, d; while (i < order / 2) { I = f(a); for (j = a, k = 0; j <= b; j = j + h, k++) x[k] = j; m = x.length; for (j = 1; j < m - 1; j++) { I += (((j % 2) !== 0) ? 4 : 2) * f(x[j]); } I = (h / 3) * (I + f(b)); g[i] = I; h /= 2; i++; } a1 = g.length; m = 1; while (a1 !== 1) { for (j = 0; j < a1 - 1; j++) h1[j] = ((Math.pow(4, m)) * g[j + 1] - g[j]) / (Math.pow(4, m) - 1); a1 = h1.length; g = h1; h1 = []; m++; } return g; }, richardson: function richardson(X, f, x, h) { function pos(X, x) { var i = 0; var n = X.length; var p; for (; i < n; i++) if (X[i] === x) p = i; return p; } var n = X.length, h_min = Math.abs(x - X[pos(X, x) + 1]), i = 0, g = [], h1 = [], y1, y2, m, a, j; while (h >= h_min) { y1 = pos(X, x + h); y2 = pos(X, x); g[i] = (f[y1] - 2 * f[y2] + f[2 * y2 - y1]) / (h * h); h /= 2; i++; } a = g.length; m = 1; while (a != 1) { for (j = 0; j < a - 1; j++) h1[j] = ((Math.pow(4, m)) * g[j + 1] - g[j]) / (Math.pow(4, m) - 1); a = h1.length; g = h1; h1 = []; m++; } return g; }, simpson: function simpson(f, a, b, n) { var h = (b - a) / n; var I = f(a); var x = []; var j = a; var k = 0; var i = 1; var m; for (; j <= b; j = j + h, k++) x[k] = j; m = x.length; for (; i < m - 1; i++) { I += ((i % 2 !== 0) ? 4 : 2) * f(x[i]); } return (h / 3) * (I + f(b)); }, hermite: function hermite(X, F, dF, value) { var n = X.length; var p = 0; var i = 0; var l = []; var dl = []; var A = []; var B = []; var j; for (; i < n; i++) { l[i] = 1; for (j = 0; j < n; j++) { if (i != j) l[i] *= (value - X[j]) / (X[i] - X[j]); } dl[i] = 0; for (j = 0; j < n; j++) { if (i != j) dl[i] += 1 / (X [i] - X[j]); } A[i] = (1 - 2 * (value - X[i]) * dl[i]) * (l[i] * l[i]); B[i] = (value - X[i]) * (l[i] * l[i]); p += (A[i] * F[i] + B[i] * dF[i]); } return p; }, lagrange: function lagrange(X, F, value) { var p = 0; var i = 0; var j, l; var n = X.length; for (; i < n; i++) { l = F[i]; for (j = 0; j < n; j++) { // calculating the lagrange polynomial L_i if (i != j) l *= (value - X[j]) / (X[i] - X[j]); } // adding the lagrange polynomials found above p += l; } return p; }, cubic_spline: function cubic_spline(X, F, value) { var n = X.length; var i = 0, j; var A = []; var B = []; var alpha = []; var c = []; var h = []; var b = []; var d = []; for (; i < n - 1; i++) h[i] = X[i + 1] - X[i]; alpha[0] = 0; for (i = 1; i < n - 1; i++) { alpha[i] = (3 / h[i]) * (F[i + 1] - F[i]) - (3 / h[i-1]) * (F[i] - F[i-1]); } for (i = 1; i < n - 1; i++) { A[i] = []; B[i] = []; A[i][i-1] = h[i-1]; A[i][i] = 2 * (h[i - 1] + h[i]); A[i][i+1] = h[i]; B[i][0] = alpha[i]; } c = jStat.multiply(jStat.inv(A), B); for (j = 0; j < n - 1; j++) { b[j] = (F[j + 1] - F[j]) / h[j] - h[j] * (c[j + 1][0] + 2 * c[j][0]) / 3; d[j] = (c[j + 1][0] - c[j][0]) / (3 * h[j]); } for (j = 0; j < n; j++) { if (X[j] > value) break; } j -= 1; return F[j] + (value - X[j]) * b[j] + jStat.sq(value-X[j]) * c[j] + (value - X[j]) * jStat.sq(value - X[j]) * d[j]; }, gauss_quadrature: function gauss_quadrature() { throw new Error('gauss_quadrature not yet implemented'); }, PCA: function PCA(X) { var m = X.length; var n = X[0].length; var flag = false; var i = 0; var j, temp1; var u = []; var D = []; var result = []; var temp2 = []; var Y = []; var Bt = []; var B = []; var C = []; var V = []; var Vt = []; for (i = 0; i < m; i++) { u[i] = jStat.sum(X[i]) / n; } for (i = 0; i < n; i++) { B[i] = []; for(j = 0; j < m; j++) { B[i][j] = X[j][i] - u[j]; } } B = jStat.transpose(B); for (i = 0; i < m; i++) { C[i] = []; for (j = 0; j < m; j++) { C[i][j] = (jStat.dot([B[i]], [B[j]])) / (n - 1); } } result = jStat.jacobi(C); V = result[0]; D = result[1]; Vt = jStat.transpose(V); for (i = 0; i < D.length; i++) { for (j = i; j < D.length; j++) { if(D[i] < D[j]) { temp1 = D[i]; D[i] = D[j]; D[j] = temp1; temp2 = Vt[i]; Vt[i] = Vt[j]; Vt[j] = temp2; } } } Bt = jStat.transpose(B); for (i = 0; i < m; i++) { Y[i] = []; for (j = 0; j < Bt.length; j++) { Y[i][j] = jStat.dot([Vt[i]], [Bt[j]]); } } return [X, D, Vt, Y]; } }); // extend jStat.fn with methods that require one argument (function(funcs) { for (var i = 0; i < funcs.length; i++) (function(passfunc) { jStat.fn[passfunc] = function(arg, func) { var tmpthis = this; // check for callback if (func) { setTimeout(function() { func.call(tmpthis, jStat.fn[passfunc].call(tmpthis, arg)); }, 15); return this; } if (typeof jStat[passfunc](this, arg) === 'number') return jStat[passfunc](this, arg); else return jStat(jStat[passfunc](this, arg)); }; }(funcs[i])); }('add divide multiply subtract dot pow exp log abs norm angle'.split(' '))); }(this.jStat, Math)); (function(jStat, Math) { var slice = [].slice; var isNumber = jStat.utils.isNumber; // flag==true denotes use of sample standard deviation // Z Statistics jStat.extend({ // 2 different parameter lists: // (value, mean, sd) // (value, array, flag) zscore: function zscore() { var args = slice.call(arguments); if (isNumber(args[1])) { return (args[0] - args[1]) / args[2]; } return (args[0] - jStat.mean(args[1])) / jStat.stdev(args[1], args[2]); }, // 3 different paramter lists: // (value, mean, sd, sides) // (zscore, sides) // (value, array, sides, flag) ztest: function ztest() { var args = slice.call(arguments); if (args.length === 4) { if(isNumber(args[1])) { var z = jStat.zscore(args[0],args[1],args[2]) return (args[3] === 1) ? (jStat.normal.cdf(-Math.abs(z),0,1)) : (jStat.normal.cdf(-Math.abs(z),0,1)* 2); } var z = args[0] return (args[2] === 1) ? (jStat.normal.cdf(-Math.abs(z),0,1)) : (jStat.normal.cdf(-Math.abs(z),0,1)*2); } var z = jStat.zscore(args[0],args[1],args[3]) return (args[1] === 1) ? (jStat.normal.cdf(-Math.abs(z), 0, 1)) : (jStat.normal.cdf(-Math.abs(z), 0, 1)*2); } }); jStat.extend(jStat.fn, { zscore: function zscore(value, flag) { return (value - this.mean()) / this.stdev(flag); }, ztest: function ztest(value, sides, flag) { var zscore = Math.abs(this.zscore(value, flag)); return (sides === 1) ? (jStat.normal.cdf(-zscore, 0, 1)) : (jStat.normal.cdf(-zscore, 0, 1) * 2); } }); // T Statistics jStat.extend({ // 2 parameter lists // (value, mean, sd, n) // (value, array) tscore: function tscore() { var args = slice.call(arguments); return (args.length === 4) ? ((args[0] - args[1]) / (args[2] / Math.sqrt(args[3]))) : ((args[0] - jStat.mean(args[1])) / (jStat.stdev(args[1], true) / Math.sqrt(args[1].length))); }, // 3 different paramter lists: // (value, mean, sd, n, sides) // (tscore, n, sides) // (value, array, sides) ttest: function ttest() { var args = slice.call(arguments); var tscore; if (args.length === 5) { tscore = Math.abs(jStat.tscore(args[0], args[1], args[2], args[3])); return (args[4] === 1) ? (jStat.studentt.cdf(-tscore, args[3]-1)) : (jStat.studentt.cdf(-tscore, args[3]-1)*2); } if (isNumber(args[1])) { tscore = Math.abs(args[0]) return (args[2] == 1) ? (jStat.studentt.cdf(-tscore, args[1]-1)) : (jStat.studentt.cdf(-tscore, args[1]-1) * 2); } tscore = Math.abs(jStat.tscore(args[0], args[1])) return (args[2] == 1) ? (jStat.studentt.cdf(-tscore, args[1].length-1)) : (jStat.studentt.cdf(-tscore, args[1].length-1) * 2); } }); jStat.extend(jStat.fn, { tscore: function tscore(value) { return (value - this.mean()) / (this.stdev(true) / Math.sqrt(this.cols())); }, ttest: function ttest(value, sides) { return (sides === 1) ? (1 - jStat.studentt.cdf(Math.abs(this.tscore(value)), this.cols()-1)) : (jStat.studentt.cdf(-Math.abs(this.tscore(value)), this.cols()-1)*2); } }); // F Statistics jStat.extend({ // Paramter list is as follows: // (array1, array2, array3, ...) // or it is an array of arrays // array of arrays conversion anovafscore: function anovafscore() { var args = slice.call(arguments), expVar, sample, sampMean, sampSampMean, tmpargs, unexpVar, i, j; if (args.length === 1) { tmpargs = new Array(args[0].length); for (i = 0; i < args[0].length; i++) { tmpargs[i] = args[0][i]; } args = tmpargs; } // 2 sample case if (args.length === 2) { return jStat.variance(args[0]) / jStat.variance(args[1]); } // Builds sample array sample = new Array(); for (i = 0; i < args.length; i++) { sample = sample.concat(args[i]); } sampMean = jStat.mean(sample); // Computes the explained variance expVar = 0; for (i = 0; i < args.length; i++) { expVar = expVar + args[i].length * Math.pow(jStat.mean(args[i]) - sampMean, 2); } expVar /= (args.length - 1); // Computes unexplained variance unexpVar = 0; for (i = 0; i < args.length; i++) { sampSampMean = jStat.mean(args[i]); for (j = 0; j < args[i].length; j++) { unexpVar += Math.pow(args[i][j] - sampSampMean, 2); } } unexpVar /= (sample.length - args.length); return expVar / unexpVar; }, // 2 different paramter setups // (array1, array2, array3, ...) // (anovafscore, df1, df2) anovaftest: function anovaftest() { var args = slice.call(arguments), df1, df2, n, i; if (isNumber(args[0])) { return 1 - jStat.centralF.cdf(args[0], args[1], args[2]); } anovafscore = jStat.anovafscore(args); df1 = args.length - 1; n = 0; for (i = 0; i < args.length; i++) { n = n + args[i].length; } df2 = n - df1 - 1; return 1 - jStat.centralF.cdf(anovafscore, df1, df2); }, ftest: function ftest(fscore, df1, df2) { return 1 - jStat.centralF.cdf(fscore, df1, df2); } }); jStat.extend(jStat.fn, { anovafscore: function anovafscore() { return jStat.anovafscore(this.toArray()); }, anovaftes: function anovaftes() { var n = 0; var i; for (i = 0; i < this.length; i++) { n = n + this[i].length; } return jStat.ftest(this.anovafscore(), this.length - 1, n - this.length); } }); // Error Bounds jStat.extend({ // 2 different parameter setups // (value, alpha, sd, n) // (value, alpha, array) normalci: function normalci() { var args = slice.call(arguments), ans = new Array(2), change; if (args.length === 4) { change = Math.abs(jStat.normal.inv(args[1] / 2, 0, 1) * args[2] / Math.sqrt(args[3])); } else { change = Math.abs(jStat.normal.inv(args[1] / 2, 0, 1) * jStat.stdev(args[2]) / Math.sqrt(args[2].length)); } ans[0] = args[0] - change; ans[1] = args[0] + change; return ans; }, // 2 different parameter setups // (value, alpha, sd, n) // (value, alpha, array) tci: function tci() { var args = slice.call(arguments), ans = new Array(2), change; if (args.length === 4) { change = Math.abs(jStat.studentt.inv(args[1] / 2, args[3] - 1) * args[2] / Math.sqrt(args[3])); } else { change = Math.abs(jStat.studentt.inv(args[1] / 2, args[2].length - 1) * jStat.stdev(args[2], true) / Math.sqrt(args[2].length)); } ans[0] = args[0] - change; ans[1] = args[0] + change; return ans; }, significant: function significant(pvalue, alpha) { return pvalue < alpha; } }); jStat.extend(jStat.fn, { normalci: function normalci(value, alpha) { return jStat.normalci(value, alpha, this.toArray()); }, tci: function tci(value, alpha) { return jStat.tci(value, alpha, this.toArray()); } }); }(this.jStat, Math));