Merge branch '1.0'

Conflicts:
	package.json

dist/jstat.js

@@ -710,6 +710,7 @@ jStat.quantiles = function quantiles(arr, quantilesArray, alphap, betap) {
   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.
@@ -733,6 +734,25 @@ jStat.percentileOfScore = function percentileOfScore(arr, score, kind) {
   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);
@@ -755,6 +775,29 @@ jStat.corrcoeff = function corrcoeff(arr1, arr2) {
       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;
 
@@ -833,7 +876,8 @@ var jProto = jStat.prototype;
     };
   })(funcs[i]);
 })(('sum sumsqrd sumsqerr product min max mean meansqerr geomean median diff ' +
-    'mode range variance stdev meandev meddev coeffvar quartiles').split(' '));
+    'mode range variance stdev meandev meddev coeffvar quartiles histogram ' +
+    'skewness kurtosis').split(' '));
 
 
 // Extend jProto with functions that take arguments. Operations on matrices are
@@ -972,7 +1016,7 @@ jStat.gammap = function gammap(a, x) {
     for (; i <= ITMAX; i++) {
       sum += del *= x / ++ap;
     }
-    return sum * Math.exp(-x + a * Math.log(x) - (aln));
+    return (sum * Math.exp(-x + a * Math.log(x) - (aln))) * jStat.gammafn(a);
   }
 
   for (; i <= ITMAX; i++) {
@@ -984,7 +1028,7 @@ jStat.gammap = function gammap(a, x) {
     h *= d * c;
   }
 
-  return 1 - h * Math.exp(-x + a * Math.log(x) - (aln));
+  return (1 - h * Math.exp(-x + a * Math.log(x) - (aln))) * jStat.gammafn(a);
 };
 
 
@@ -1807,7 +1851,7 @@ jStat.extend(jStat.normal, {
 // extend pareto function with static methods
 jStat.extend(jStat.pareto, {
   pdf: function pdf(x, scale, shape) {
-    if (x <= scale)
+    if (x < scale)
       return undefined;
     return (shape * Math.pow(scale, shape)) / Math.pow(x, shape + 1);
   },

doc/md/vector.md

@@ -606,10 +606,26 @@ If pass boolean true as first argument, then return median absolute deviation of
 
 And the two can be combined.
 
-    jStat[[1,2],[3,5]]).meddev(true,function( result ) {
+    jStat([[1,2],[3,5]]).meddev(true,function( result ) {
         // result === 0.25
     });
 
+### skewness()
+
+**skewness( array )**
+
+Return the skewness of a vector (third standardized moment).
+
+    jStat.skewness([1,2,2,3,5]) === 0.75003...
+
+### kurtosis()
+
+**kurtosis( array )**
+
+Return the excess kurtosis of a vector (fourth standardized moment - 3).
+
+    jStat.kurtosis([1,2,3,4]) === -0.63610...
+
 ### coeffvar()
 
 **coeffvar( array )**
@@ -688,6 +704,16 @@ less or equal than (if `kind == 'weak'`) score. Default is `'weak'`.
 
      jStat.percentileOfScore([1, 2, 3, 4, 5, 6], 3), 0.5, 'weak') === 0.5;
 
+### histogram()
+
+**histogram( dataArray, [numBins] )**
+
+The histogram data defined as the number of dataArray elements found in 
+equally sized bins across the range of dataArray. Default number
+of bins is 4.
+
+     jStat.histogram([100, 101, 102, 230, 304, 305, 400], 3) === [3, 1, 3];
+
 ### covariance()
 
 **covariance( array, array )**

src/distribution.js

@@ -128,6 +128,12 @@ jStat.extend(jStat.centralF, {
       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));
@@ -462,7 +468,7 @@ jStat.extend(jStat.normal, {
 // extend pareto function with static methods
 jStat.extend(jStat.pareto, {
   pdf: function pdf(x, scale, shape) {
-    if (x <= scale)
+    if (x < scale)
       return undefined;
     return (shape * Math.pow(scale, shape)) / Math.pow(x, shape + 1);
   },

src/special.js

@@ -89,7 +89,7 @@ jStat.gammap = function gammap(a, x) {
     for (; i <= ITMAX; i++) {
       sum += del *= x / ++ap;
     }
-    return sum * Math.exp(-x + a * Math.log(x) - (aln));
+    return (sum * Math.exp(-x + a * Math.log(x) - (aln))) * jStat.gammafn(a);
   }
 
   for (; i <= ITMAX; i++) {
@@ -101,7 +101,7 @@ jStat.gammap = function gammap(a, x) {
     h *= d * c;
   }
 
-  return 1 - h * Math.exp(-x + a * Math.log(x) - (aln));
+  return (1 - h * Math.exp(-x + a * Math.log(x) - (aln))) * jStat.gammafn(a);
 };
 
 

src/vector.js

@@ -248,6 +248,7 @@ jStat.quantiles = function quantiles(arr, quantilesArray, alphap, betap) {
   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.
@@ -271,6 +272,25 @@ jStat.percentileOfScore = function percentileOfScore(arr, score, kind) {
   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);
@@ -293,6 +313,29 @@ jStat.corrcoeff = function corrcoeff(arr1, arr2) {
       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;
 
@@ -371,7 +414,8 @@ var jProto = jStat.prototype;
     };
   })(funcs[i]);
 })(('sum sumsqrd sumsqerr product min max mean meansqerr geomean median diff ' +
-    'mode range variance stdev meandev meddev coeffvar quartiles').split(' '));
+    'mode range variance stdev meandev meddev coeffvar quartiles histogram ' +
+    'skewness kurtosis').split(' '));
 
 
 // Extend jProto with functions that take arguments. Operations on matrices are

test/distribution/beta-test.js

@@ -0,0 +1,45 @@
+var vows = require('vows');
+var assert = require('assert');
+var suite = vows.describe('jStat.distribution');
+
+require('../env.js');
+
+suite.addBatch({
+  'beta pdf': {
+    'topic': function() {
+      return jStat;
+    },
+    'check pdf calculation': function(jStat) {
+      // Non-log form of the Beta pdf
+      function pdf(x, alpha, beta) {
+        if (x > 1 || x < 0)
+          return 0;
+        return (Math.pow(x, alpha - 1) * Math.pow(1 - x, beta - 1)) /
+            jStat.betafn(alpha, beta);
+      }
+
+      var tol = 0.0000001;
+      var args = [0, 0.1, 0.3, 0.5, 0.7, 0.9, 1];
+      var arg;
+
+      for (var i = 0; i < args.length; i++) {
+        arg = args[i];
+        assert.epsilon(tol, jStat.beta.pdf(arg, 0.1, 0.1), pdf(arg, 0.1, 0.1));
+        assert.epsilon(tol, jStat.beta.pdf(arg, 1, 1), pdf(arg, 1, 1));
+        assert.epsilon(tol, jStat.beta.pdf(arg, 10, 50), pdf(arg, 10, 50));
+
+        // Show that the log form of the pdf performs better for
+        // large parameter values.
+        assert(!isNaN(jStat.beta.pdf(arg, 1000, 5000)),
+               'New Beta pdf is NaN for large parameter values.');
+        assert(isNaN(pdf(arg, 1000, 5000)),
+               'Old Beta pdf is not NaN for large parameter values.');
+      }
+
+      assert.equal(jStat.beta.pdf(0, 1, 4), 4);
+      assert.equal(jStat.beta.pdf(1, 4, 1), 4);
+    }
+  }
+});
+
+suite.export(module);

test/distribution/binomial-test.js

@@ -0,0 +1,27 @@
+var vows = require('vows');
+var assert = require('assert');
+var suite = vows.describe('jStat.distribution');
+
+require('../env.js');
+
+suite.addBatch({
+  'binomial pdf': {
+    'topic': function() {
+      return jStat;
+    },
+    //checked against R's dbinom(k, n, p)
+    'check pdf calculation': function(jStat) {
+      var tol = 0.0000001;
+      assert.epsilon(tol, jStat.binomial.pdf(10, 25, 0.5), 0.09741664);
+      assert.epsilon(tol, jStat.binomial.pdf(50, 1000, 0.05), 0.05778798);
+    },
+    //Checked against r's pbinom(k, n, p)
+    'check cdf calculation': function(jStat) {
+      var tol = 0.0000001;
+      assert.epsilon(tol, jStat.binomial.cdf(10, 25, 0.5), 0.2121781);
+      assert.epsilon(tol, jStat.binomial.cdf(50, 1000, 0.05), 0.537529);
+    }
+  }
+});
+
+suite.export(module);

test/distribution/central-f-test.js

@@ -0,0 +1,44 @@
+var vows = require('vows');
+var assert = require('assert');
+var suite = vows.describe('jStat.distribution');
+
+require('../env.js');
+
+suite.addBatch({
+  'central F distribution': {
+    'topic': function() {
+      return jStat;
+    },
+    //Check against R's df(x, df1, df2)
+    'check pdf calculation': function(jStat) {
+      var tol = 0.0000001;
+
+      var zeroth = jStat.centralF.pdf(0.2, 1, 3);
+      assert.epsilon(tol, zeroth, 0.722349);
+
+      var first = jStat.centralF.pdf(1, 100, 100);
+      assert.epsilon(tol, first, 1.989731);
+
+      var second = jStat.centralF.pdf(2.5, 50, 200);
+      assert.epsilon(tol, second, 0.00003610325);
+
+      var third = jStat.centralF.pdf(0.8, 2, 10);
+      assert.epsilon(tol, third, 0.4104423);
+
+      var fourth = jStat.centralF.pdf(0.4, 3, 10);
+      assert.epsilon(tol, fourth, 0.6733766);
+
+      var first_at_zero = jStat.centralF.pdf(0.0, 3, 5);
+      assert.epsilon(tol, first_at_zero, 0);
+
+      var second_at_zero = jStat.centralF.pdf(0.0, 2, 1);
+      assert.epsilon(tol, second_at_zero, 1);
+
+      var third_at_zero = jStat.centralF.pdf(0.0, 1, 1);
+      assert.strictEqual(third_at_zero, Infinity);
+
+    }
+  },
+});
+
+suite.export(module);