Repository: josdejong/mathjs
Pull request #366 (2015-05-12 10:08:52)

Basic quantile implementation and tests (see issues #305 and #354)
There is also an implementation that assumes the input is sorted in ascending order, to maintain an O(1) runtime.


Patch:
diff --git a/docs/expressions/customization.md b/docs/expressions/customization.md
index 6dc28193a..c680e94a5 100644
--- a/docs/expressions/customization.md
+++ b/docs/expressions/customization.md
@@ -10,7 +10,7 @@ return value by writing a *transform* for the function. A transform is a
 function wrapping around a function to be transformed or completely replaces
 a function.
 
-For example, the functions for math.js use zero-based matrix indices (as is
+For example, the functions or math.js use zero-based matrix indices (as is
 common in programing languages), but the expression parser uses one-based
 indices. To enable this, all functions dealing with indices have a transform,
 which changes input from one-based to zero-based, and transforms output (and
diff --git a/lib/function/statistics/quantileSeq.js b/lib/function/statistics/quantileSeq.js
new file mode 100644
index 000000000..15e8d6f52
--- /dev/null
+++ b/lib/function/statistics/quantileSeq.js
@@ -0,0 +1,227 @@
+'use strict';
+
+module.exports = function (math) {
+  var util = require('../../util/index');
+  var collection = math.collection;
+  var BigNumber = math.type.BigNumber;
+  var Unit = require('../../type/Unit');
+
+  var isArray = Array.isArray;
+  var isBoolean = util.boolean.isBoolean;
+  var isCollection = collection.isCollection;
+  var isInteger = util.number.isInteger;
+  var isNumber = util.number.isNumber;
+
+  var add = math.add;
+  var flatten = util.array.flatten;
+  var multiply = math.multiply;
+
+  /**
+   * Compute the prob order quantile of a matrix or a list with values.
+   * The sequence is sorted and the middle value is returned.
+   * Supported types of sequence values are: Number, BigNumber, Unit
+   * Supported types of probablity are: Number, BigNumber
+   *
+   * In case of a (multi dimensional) array or matrix, the prob order quantile
+   * of all elements will be calculated.
+   *
+   * Syntax:
+   *
+   *     math.quantileSeq(A, prob[, sorted])
+   *     math.quantileSeq(A, [prob1, prob2, ...][, sorted])
+   *     math.quantileSeq(A, N[, sorted])
+   *
+   * Examples:
+   *
+   *     math.quantileSeq([3, -1, 5, 7], 0.5);         // returns 4
+   *     math.quantileSeq([3, -1, 5, 7], [1/3, 2/3]);  // returns [3, 5]
+   *     math.quantileSeq([3, -1, 5, 7], 2);           // returns [3, 5]
+   *     math.quantileSeq([-1, 3, 5, 7], 0.5, true);   // returns 4
+   *
+   * See also:
+   *
+   *     median, mean, min, max, sum, prod, std, var
+   *
+   * @param {Array, Matrix} data                A single matrix or Array
+   * @param {Number, BigNumber, Array} probOrN  prob is the order of the quantile, while N is
+   *                                            the amount of evenly distributed steps of
+   *                                            probabilities; only one of these options can
+   *                                            be provided
+   * @param {Boolean} sorted=False              is data sorted in ascending order
+   * @return {Number, BigNumber, Unit, Array}   Quantile(s)
+   */
+  math.quantileSeq = function quantileSeq(data, probOrN, sorted) {
+    var probArr, dataArr, one;
+
+    if (arguments.length < 2 || arguments.length > 3) {
+      throw new SyntaxError('Function quantileSeq requires two or three parameters');
+    }
+
+    if (isCollection(data)) {
+      sorted = sorted || false;
+      if (isBoolean(sorted)) {
+        dataArr = data.valueOf();
+        if (isNumber(probOrN)) {
+          if (probOrN < 0) {
+            throw new Error('N/prob must be non-negative');
+          }
+
+          if (probOrN <= 1) {
+            // quantileSeq([a, b, c, d, ...], prob[,sorted])
+            return _quantileSeq(dataArr, probOrN, sorted);
+          }
+
+          if (probOrN > 1) {
+            // quantileSeq([a, b, c, d, ...], N[,sorted])
+            if (!isInteger(probOrN)) {
+              throw new Error('N must be a positive integer');
+            }
+
+            var nPlusOne = probOrN + 1;
+            probArr = new Array(probOrN);
+            for (var i = 0; i < probOrN;) {
+              probArr[i] = _quantileSeq(dataArr, (++i) / nPlusOne, sorted);
+            }
+            return probArr;
+          }
+        }
+
+        if (probOrN instanceof BigNumber) {
+          if (probOrN.isNegative()) {
+            throw new Error('N/prob must be non-negative');
+          }
+
+          one = probOrN.constructor.ONE;
+
+          if (probOrN.lte(one)) {
+            // quantileSeq([a, b, c, d, ...], prob[,sorted])
+            return _quantileSeq(dataArr, probOrN, sorted);
+          }
+
+          if (probOrN.gt(one)) {
+            // quantileSeq([a, b, c, d, ...], N[,sorted])
+            if (!probOrN.isInteger()) {
+              throw new Error('N must be a positive integer');
+            }
+
+            var intN = probOrN.toNumber();
+            var nPlusOne = new BigNumber(intN + 1);
+
+            probArr = new Array(intN);
+            for (var i = 0; i < intN;) {
+              probArr[i] = _quantileSeq(dataArr, new BigNumber(++i).div(nPlusOne), sorted);
+            }
+            return probArr;
+          }
+        }
+
+        if (isArray(probOrN)) {
+          // quantileSeq([a, b, c, d, ...], [prob][,sorted])
+          probArr = new Array(probOrN.length);
+          for (var i = 0; i < probArr.length; ++i) {
+            var currProb = probOrN[i];
+            if (isNumber(currProb)) {
+              if (currProb < 0 || currProb > 1) {
+                throw new Error('Probability must be between 0 and 1, inclusive');
+              }
+            } else if (currProb instanceof BigNumber) {
+              one = currProb.constructor.ONE;
+              if (currProb.isNegative() || currProb.gt(one)) {
+                throw new Error('Probability must be between 0 and 1, inclusive');
+              }
+            } else {
+              throw new math.error.UnsupportedTypeError('quantileSeq', math['typeof'](currProb));
+            }
+
+            probArr[i] = _quantileSeq(dataArr, currProb, sorted);
+          }
+          return probArr;
+        }
+
+        throw new math.error.UnsupportedTypeError('quantileSeq', math['typeof'](probOrN));
+      }
+
+      throw new math.error.UnsupportedTypeError('quantileSeq', math['typeof'](sorted));
+    }
+
+    throw new math.error.UnsupportedTypeError('quantileSeq', math['typeof'](data));
+  };
+
+  /**
+   * Calculate the prob order quantile of an n-dimensional array.
+   * 
+   * @param {Array} array
+   * @param {Number, BigNumber} prob
+   * @param {Boolean} sorted
+   * @return {Number, BigNumber, Unit} prob order quantile
+   * @private
+   */
+  function _quantileSeq(array, prob, sorted) {
+    var flat = flatten(array);
+    var len = flat.length;
+    if (len === 0) {
+      throw new Error('Cannot calculate quantile of an empty sequence');
+    }
+
+    if (!sorted) {
+      flat.sort(math.compare);
+    }
+
+    if (isNumber(prob)) {
+      var index = prob * (len-1);
+      var fracPart = index % 1;
+      if (fracPart === 0) {
+        var value = flat[index];
+
+        typecheck(value);
+
+        return value;
+      }
+
+      var integerPart = Math.floor(index);
+      var left = flat[integerPart];
+      var right = flat[integerPart+1];
+
+      typecheck(left, right);
+
+      return add(multiply(left, 1 - fracPart), multiply(right, fracPart));
+    }
+
+    // If prob is a BigNumber
+    var index = prob.times(len-1);
+    if (index.isInteger()) {
+      var value = flat[index.toNumber()];
+
+      typecheck(value);
+
+      return value;
+    }
+
+    var integerPart = index.floor();
+    var fracPart = index.minus(integerPart);
+
+    var integerPartNumber = integerPart.toNumber();
+    var left = flat[integerPartNumber];
+    var right = flat[integerPartNumber+1];
+
+    typecheck(left, right);
+
+    var one = fracPart.constructor.ONE;
+    return add(multiply(left, one.minus(fracPart)), multiply(right, fracPart));
+  }
+
+  /**
+   * Check if array value types are valid, throw error otherwise.
+   * @param {*} x
+   * @param {*} y
+   * @private
+   */
+  function typecheck(x, y) {
+    if (!isNumber(x) && !(x instanceof BigNumber) && !(x instanceof Unit)) {
+      throw new math.error.UnsupportedTypeError('quantileSeq', math['typeof'](x));
+    }
+    if (y !== undefined && !isNumber(y) && !(y instanceof BigNumber) && !(y instanceof Unit)) {
+      throw new math.error.UnsupportedTypeError('quantileSeq', math['typeof'](y));
+    }
+  }
+};
diff --git a/lib/math.js b/lib/math.js
index 4dc8e7a15..c937a9848 100644
--- a/lib/math.js
+++ b/lib/math.js
@@ -153,7 +153,7 @@ function create (config) {
   if (typeof BigNumber.convert !== 'function') {
     /**
      * Try to convert a Number in to a BigNumber.
-     * If the number has 15 or mor significant digits, the Number cannot be
+     * If the number has 15 or more significant digits, the Number cannot be
      * converted to BigNumber and will return the original number.
      * @param {Number} number
      * @return {BigNumber | Number} bignumber
@@ -329,6 +329,7 @@ function create (config) {
   require('./function/statistics/mean')(math, _config);
   require('./function/statistics/median')(math, _config);
   require('./function/statistics/prod')(math, _config);
+  require('./function/statistics/quantileSeq')(math, _config);
   require('./function/statistics/std')(math, _config);
   require('./function/statistics/sum')(math, _config);
   require('./function/statistics/var')(math, _config);


Test patch:
diff --git a/test/function/statistics/quantileSeq.test.js b/test/function/statistics/quantileSeq.test.js
new file mode 100644
index 000000000..2532a4434
--- /dev/null
+++ b/test/function/statistics/quantileSeq.test.js
@@ -0,0 +1,186 @@
+var assert = require('assert'),
+    approx = require('../../../tools/approx'),
+    math = require('../../../index'),
+    bignumber = math.bignumber,
+    quantileSeq = math.quantileSeq;
+
+describe('quantileSeq', function() {
+
+  it('should return the quantileSeq from an array with number probability', function() {
+    var lst = [3.7, 2.7, 3.3, 1.3, 2.2, 3.1];
+    assert.equal(quantileSeq(lst, 0), 1.3);
+    assert.equal(quantileSeq(lst, 0.1), 1.75);
+    assert.equal(quantileSeq(lst, 0.2), 2.2);
+    assert.equal(quantileSeq(lst, 0.25), 2.325);
+    assert.equal(quantileSeq(lst, 0.25, false), 2.325);
+    assert.equal(quantileSeq(lst, 0.3), 2.45);
+    assert.equal(quantileSeq(lst, 0.4), 2.7);
+    approx.equal(quantileSeq(lst, 0.5), 2.9);
+    assert.equal(quantileSeq(lst, 0.6), 3.1);
+    assert.equal(quantileSeq(lst, 0.7), 3.2);
+    approx.equal(quantileSeq(lst, 0.75), 3.25);
+    assert.equal(quantileSeq(lst, 0.8), 3.3);
+    assert.equal(quantileSeq(lst, 0.9), 3.5);
+    assert.equal(quantileSeq(lst, 1), 3.7);
+  });
+
+  it('should return the quantileSeq from an ascending array with number probability', function() {
+    var lst = [1.3, 2.2, 2.7, 3.1, 3.3, 3.7];
+    assert.equal(quantileSeq(lst, 0, true), 1.3);
+    assert.equal(quantileSeq(lst, 0.1, true), 1.75);
+    assert.equal(quantileSeq(lst, 0.2, true), 2.2);
+    assert.equal(quantileSeq(lst, 0.25, true), 2.325);
+    assert.equal(quantileSeq(lst, 0.3, true), 2.45);
+    assert.equal(quantileSeq(lst, 0.4, true), 2.7);
+    approx.equal(quantileSeq(lst, 0.5, true), 2.9);
+    assert.equal(quantileSeq(lst, 0.6, true), 3.1);
+    assert.equal(quantileSeq(lst, 0.7, true), 3.2);
+    approx.equal(quantileSeq(lst, 0.75, true), 3.25);
+    assert.equal(quantileSeq(lst, 0.8, true), 3.3);
+    assert.equal(quantileSeq(lst, 0.9, true), 3.5);
+    assert.equal(quantileSeq(lst, 1, true), 3.7);
+  });
+
+  it('should return the quantileSeq from an array with BigNumber probability', function() {
+    var lst = [3.7, 2.7, 3.3, 1.3, 2.2, 3.1];
+    assert.equal(quantileSeq(lst, bignumber(0)), 1.3);
+    assert.equal(quantileSeq(lst, bignumber(0.1)), 1.75);
+    assert.equal(quantileSeq(lst, bignumber(0.2)), 2.2);
+    assert.equal(quantileSeq(lst, bignumber(0.25)), 2.325);
+    assert.equal(quantileSeq(lst, bignumber(0.3)), 2.45);
+    assert.equal(quantileSeq(lst, bignumber(0.4)), 2.7);
+    assert.equal(quantileSeq(lst, bignumber(0.5)), 2.9);
+    assert.equal(quantileSeq(lst, bignumber(0.6)), 3.1);
+    assert.equal(quantileSeq(lst, bignumber(0.7)), 3.2);
+    assert.equal(quantileSeq(lst, bignumber(0.75)), 3.25);
+    assert.equal(quantileSeq(lst, bignumber(0.8)), 3.3);
+    assert.equal(quantileSeq(lst, bignumber(0.9)), 3.5);
+    assert.equal(quantileSeq(lst, bignumber(1)), 3.7);
+  });
+
+  it('should return the quantileSeq of an array of bignumbers with number probability', function() {
+    approx.equal(quantileSeq([bignumber(0.5377),bignumber(1.8339),bignumber(-2.2588),bignumber(0.8622),
+                              bignumber(0.3188),bignumber(-1.3077),bignumber(-0.4336),bignumber(0.3426),
+                              bignumber(3.5784),bignumber(2.7694)],0.3),
+                              0.09308);
+  });
+
+  it('should return the quantileSeq of an array of bignumbers with BigNumber probability', function() {
+    assert.deepEqual(quantileSeq([bignumber(0.5377),bignumber(1.8339),bignumber(-2.2588),bignumber(0.8622),
+                                  bignumber(0.3188),bignumber(-1.3077),bignumber(-0.4336),bignumber(0.3426),
+                                  bignumber(3.5784),bignumber(2.7694)],bignumber(0.3)),
+                                  bignumber(0.09308));
+  });
+
+  it('should return the quantileSeq of units', function() {
+    assert.deepEqual(quantileSeq([math.unit('5mm'), math.unit('15mm'), math.unit('10mm')], 0.5), math.unit('10mm'));
+  });
+
+  it('should return the quantileSeq from an 1d matrix', function() {
+    assert.equal(quantileSeq(math.matrix([2,4,6,8,10,12,14]), 0.25), 5);
+  });
+
+  it('should return the quantileSeq from a 2d array', function() {
+    approx.equal(quantileSeq([
+      [3.7, 2.7, 3.3],
+      [1.3, 2.2, 3.1]
+    ], 0.75), 3.25);
+  });
+
+  it('should return the quantileSeq from an ascending 2d array', function() {
+    approx.equal(quantileSeq([
+      [1.3, 2.2, 2.7],
+      [3.1, 3.3, 3.7]
+    ], 0.75, true), 3.25);
+  });
+
+  it('should return the quantileSeq from a 2d matrix', function() {
+    approx.equal(quantileSeq(math.matrix([
+      [3.7, 2.7, 3.3],
+      [1.3, 2.2, 3.1]
+    ]), 0.75), 3.25);
+  });
+
+  it('should return the quantileSeq from an ascending 2d matrix', function() {
+    approx.equal(quantileSeq(math.matrix([
+      [1.3, 2.2, 2.7],
+      [3.1, 3.3, 3.7]
+    ]), 0.75, true), 3.25);
+  });
+
+  it('should return list quantiles for list of number probabilities', function() {
+    var lst = [3.7, 2.7, 3.3, 1.3, 2.2, 3.1];
+    approx.deepEqual(quantileSeq(lst, [0.25, 0.5, 0.75]), [2.325, 2.9, 3.25]);
+    approx.deepEqual(quantileSeq(lst, [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]),
+                     [1.75, 2.2, 2.45, 2.7, 2.9, 3.1, 3.2, 3.3, 3.5]);
+  });
+
+  it('should return list quantiles for list of BigNumber probabilities', function() {
+    var lst = [3.7, 2.7, 3.3, 1.3, 2.2, 3.1];
+    assert.deepEqual(quantileSeq(lst, [bignumber(0.25), bignumber(0.5), bignumber(0.75)]),
+                     [bignumber(2.325), bignumber(2.9), bignumber(3.25)]);
+    assert.equal(quantileSeq(lst, [bignumber(0.1), bignumber(0.2), bignumber(0.3), bignumber(0.4),
+                                   bignumber(0.5), bignumber(0.6), bignumber(0.7), bignumber(0.8),
+                                   bignumber(0.9)]).toString(), '1.75,2.2,2.45,2.7,2.9,3.1,3.2,3.3,3.5');
+  });
+
+  it('should return list quantiles for list of number and BigNumber probabilities', function() {
+    var lst = [3.7, 2.7, 3.3, 1.3, 2.2, 3.1];
+    approx.deepEqual(quantileSeq(lst, [0.25, bignumber(0.5), 0.75]), [2.325, 2.9, 3.25]);
+    approx.deepEqual(quantileSeq(lst, [0.1, 0.2, bignumber(0.3), 0.4, 0.5, 0.6, 0.7, bignumber(0.8), 0.9]),
+                     [1.75, 2.2, 2.45, 2.7, 2.9, 3.1, 3.2, 3.3, 3.5]);
+  });
+
+  it('should return the evenly number spaced quantiles of an array', function() {
+    var lst = [3.7, 2.7, 3.3, 1.3, 2.2, 3.1];
+    approx.deepEqual(quantileSeq(lst, 3), [2.325, 2.9, 3.25]);
+    approx.deepEqual(quantileSeq(lst, 9), [1.75, 2.2, 2.45, 2.7, 2.9, 3.1, 3.2, 3.3, 3.5]);
+  });
+
+  it('should return the evenly BigNumber spaced quantiles of an array', function() {
+    var lst = [3.7, 2.7, 3.3, 1.3, 2.2, 3.1];
+    assert.deepEqual(quantileSeq(lst, bignumber(3)), [bignumber(2.325), bignumber(2.9), bignumber(3.25)]);
+    assert.equal(quantileSeq(lst, bignumber(9)).toString(), "1.75,2.2,2.45,2.7,2.9,3.1,3.2,3.3,3.5");
+  });
+
+  it('should throw an error if called with invalid number of arguments', function() {
+    assert.throws(function() {quantileSeq()}, SyntaxError);
+    assert.throws(function() {quantileSeq(2)}, SyntaxError);
+    assert.throws(function() {quantileSeq([], 2, 3, 1)}, SyntaxError);
+  });
+
+  it('should throw an error if called with unsupported type of arguments', function() {
+    assert.throws(function () {quantileSeq([2,4,6,8,10,12,14], 0.25, 10)}, math.error.UnsupportedTypeError);
+    assert.throws(function () {quantileSeq([2,4,6,8,10,12,14], [0.25, 2])}, math.error.UnsuppoError);
+    assert.throws(function () {quantileSeq('A', 'C', 'B')}, math.error.UnsupportedTypeError);
+    assert.throws(function () {quantileSeq(true, false, true)}, math.error.UnsupportedTypeError);
+    assert.throws(function () {quantileSeq(0, 'B')}, math.error.UnsupportedTypeError);
+    assert.throws(function () {quantileSeq(math.complex(2,3), math.complex(-1,2))}, TypeError);
+  });
+
+  it('should throw error for bad probabilities and splits', function() {
+    assert.throws(function () {quantileSeq([2,4,6,8,10,12,14], [0.23, 2, 0.2])}, Error);
+    assert.throws(function () {quantileSeq([2,4,6,8,10,12,14], [0.23, bignumber(2), 0.2])}, Error);
+    assert.throws(function () {quantileSeq([2,4,6,8,10,12,14], -2)}, Error);
+    assert.throws(function () {quantileSeq([2,4,6,8,10,12,14], bignumber(-2))}, Error);
+  });
+
+  it('should throw an error if called with an empty array', function() {
+    assert.throws(function() {quantileSeq([])});
+  });
+
+  it('should not mutate the input', function () {
+    var a = [3,2,1];
+    quantileSeq(a, 0.2);
+    quantileSeq(a, 2);
+    quantileSeq(a, [0.1, 0.3]);
+    assert.deepEqual(a,[3,2,1]);
+  });
+
+  /*
+  it('should LaTeX quantileSeq', function () {
+    var expression = math.parse('quantileSeq(1,2,3,4,0.3)');
+    assert.equal(expression.toTex(), '\\mathrm{quantile}\\left(1,2,3,4,0.3\\right)');
+  });
+  */
+});