# thumbtack/abba

overhaul stats code, improve hypothesis test and use jStat

```- use jStat's normal pdf/cdf/inverse cdf functions
- new implementation of the p-value computation that is better and
more statistically valid```
1 parent 689d568 commit 3194e6d5fd14d4f827cd135aec9610b00da747cd showard committed Feb 24, 2012
Showing with 291 additions and 129 deletions.
1. +170 −84 abtest/stats.js
2. +114 −45 abtest/tests/stats_test.js
3. +1 −0 abtest/tests/test_runner.html
4. +1 −0 demo/abtest.html
5. +5 −0 lib/jstat-min-08fb5e8.js
254 abtest/stats.js
 @@ -1,78 +1,92 @@ var ABTest = {}; -/* Polynomial and rational approximations to standard normal probability functions. From: - - Abramowitz, Milton; Stegun, Irene A., eds. (1972), Handbook of Mathematical Functions with - Formulas, Graphs, and Mathematical Tables, New York: Dover Publications, ISBN 978-0-486-61272-0 - - Available online at http://people.math.sfu.ca/~cbm/aands/ -*/ -ABTest.NormalDistribution = function() {} +// Friendly wrapper over jStat's normal distribution functions. +ABTest.NormalDistribution = function(mean, standardDeviation) { + if (mean === undefined) { + mean = 0; + } + if (standardDeviation === undefined) { + standardDeviation = 1; + } + this.mean = mean; + this.standardDeviation = standardDeviation; +}; ABTest.NormalDistribution.prototype = { - density: function(zValue) { - return 1 / Math.sqrt(2 * Math.PI) * Math.exp(-zValue * zValue / 2); - }, - - // Returns P(x < zValue) for x standard normal. zValue may be any number. - cdf: function(zValue) { - // Formula 26.2.17, http://people.math.sfu.ca/~cbm/aands/page932.htm - // Valid for zValue >= 0, abs(error) < 7.5 x 10^-8 - var p = 0.2316419; - var b1 = 0.319381530; - var b2 = -0.356563782; - var b3 = 1.781477937; - var b4 = -1.821255978; - var b5 = 1.330274429; - - var isInverted = false; - if (zValue < 0) { - zValue = -zValue; - isInverted = true; - } + density: function(value) { + return jStat.normal.pdf(value, this.mean, this.standardDeviation); + }, - var t = 1 / (1 + p * zValue); - var density = this.density(zValue); - var probability = 1 - density * t * (b1 + t * (b2 + t * (b3 + t * (b4 + t * b5)))); - if (isInverted) { - probability = 1 - probability; - } - return probability; + // Returns P(x < value) for x standard normal. value may be any number. + cdf: function(value) { + return jStat.normal.cdf(value, this.mean, this.standardDeviation); + }, + + // Returns P(x > value) for x standard normal. value may be any number. + survival: function(value) { + return 1 - this.cdf(value); }, - // Returns P(x > zValue) for x standard normal. zValue may be any number. - survival: function(zValue) { - return 1 - this.cdf(zValue); + // Returns z such that P(x < z) = probability for x standard normal. + // probability must be in (0, 1). + inverseCdf: function(probability) { + return jStat.normal.inv(probability, this.mean, this.standardDeviation); }, // Returns z such that P(x > z) = probability for x standard normal. // probability must be in (0, 1). inverseSurvival: function(probability) { - // Formula 26.2.23, http://people.math.sfu.ca/~cbm/aands/page933.htm - // Valid for 0 < probability <= 0.5, abs(error) < 4.5 x 10^-4 - var c0 = 2.515517; - var c1 = 0.802853; - var c2 = 0.010328; - var d1 = 1.432788; - var d2 = 0.189269; - var d3 = 0.001308; - - var multiplier = 1; - if (probability > 0.5) { - probability = 1 - probability; - multiplier = -1; + return this.mean - (this.inverseCdf(probability) - this.mean); + }, +}; + +/* Distribution functions for the binomial distribution. Relies entirely on the normal + approximation. + + jStat's binomial functions do not seem be to reliable or as performant for large cases. This + class could be improved by making it compute exact binomial functions for small cases and fall + back to the normal approximation for large cases. +*/ +ABTest.BinomialDistribution = function(numSamples, probability) { + this.numSamples = numSamples; + this.probability = probability; + this.expectation = numSamples * probability; + this.standardDeviation = Math.sqrt(this.expectation * (1 - probability)); + + // normal approximation to this binomial distribution + this._normal = new ABTest.NormalDistribution(this.expectation, this.standardDeviation); + this._lowerTailProbability = this._normal.cdf(-0.5); + this._upperTailProbability = this._normal.survival(numSamples + 0.5); +}; +ABTest.BinomialDistribution.prototype = { + mass: function(count) { + return this._normal.density(count); + }, + + _rescaleProbability: function(probability) { + return probability / (1 - this._lowerTailProbability - this._upperTailProbability); + }, + + cdf: function(count) { + if (count < 0) { + return 0; + } else if (count >= this.numSamples) { + return 1; + } else { + return this._rescaleProbability( + this._normal.cdf(count + 0.5) - this._lowerTailProbability); } + }, - var t = Math.sqrt( - Math.log(1 / (probability * probability)) - ); - var zEstimate = t - (c0 + t * (c1 + t * c2)) / (1 + t * (d1 + t * (d2 + t * d3))); - return zEstimate * multiplier; + survival: function(count) { + return 1 - this.cdf(count); }, - // Returns z such that P(x < z) = probability for x standard normal. - // probability must be in (0, 1). inverseCdf: function(probability) { - return -this.inverseSurvival(probability); + return Math.max(0, Math.min(this.numSamples, this._normal.inverseCdf(probability))); + }, + + inverseSurvival: function(probability) { + return Math.max(0, Math.min(this.numSamples, this._normal.inverseSurvival(probability))); }, }; @@ -132,7 +146,7 @@ ABTest.Proportion.prototype = { ABTest.ProportionComparison = function(baseline, trial) { this.baseline = baseline; this.trial = trial; - this._normal = new ABTest.NormalDistribution(); + this._standardNormal = new ABTest.NormalDistribution(); } ABTest.ProportionComparison.prototype = { // Generate an estimate of the difference in success rates between the trial and the baseline. @@ -152,40 +166,112 @@ ABTest.ProportionComparison.prototype = { return new ABTest.ValueWithError(ratio, error); }, - /* Perform a large-sample z-test of null hypothesis H0: pBaseline == pTrial against - alternative hypothesis H1: pBaseline < pTrial. Return the (one-tailed) p-value. + /* Compute various values useful for comparing proportions with null hypothesis that they have + the same probability of success + */ + _computeTestValues: function() { + var pooledProportion = + (this.baseline.numSuccesses + this.trial.numSuccesses) + / (this.baseline.numSamples + this.trial.numSamples); + var expectedDifference = + pooledProportion * (this.trial.numSamples - this.baseline.numSamples); + var observedDifference = this.trial.numSuccesses - this.baseline.numSuccesses; + return { + pooledProportion: pooledProportion, + expectedDifference: expectedDifference, + varianceOfDifference: + pooledProportion * (1 - pooledProportion) + * (this.baseline.numSamples + this.trial.numSamples), + observedAbsoluteDeviation: Math.abs(observedDifference - expectedDifference), + }; + }, + + /* For the given binomial distribution, compute an interval that covers at least + (1 - coverageAlpha) of the total probability mass, centered at the expectation (unless we're + at the boundary). Uses the normal approximation. + */ + _binomialCoverageInterval: function(distribution, coverageAlpha) { + if (distribution.numSamples < 1000) { + // don't even bother trying to optimize for small-ish sample sizes + return [0, distribution.numSamples]; + } else { + return [ + Math.floor(distribution.inverseCdf(coverageAlpha / 2)), + Math.ceil(distribution.inverseSurvival(coverageAlpha / 2)), + ]; + } + }, + + /* Given the probability of an event, compute the probability that it happens at least once in + numTimes independent trials. This is used to adjust a p-value for multiple comparisons. + When used to adjust alpha instead, this is called a Sidak correction (the logic is the same, + the formula is inverted): + + http://en.wikipedia.org/wiki/Bonferroni_correction#.C5.A0id.C3.A1k_correction + */ + _probabilityUnion: function(probability, numTimes) { + return 1 - Math.pow(1 - probability, numTimes); + }, - zMultiplier: test z-value will be multiplied by this factor before computing a p-value. + /* Compute a p-value testing null hypothesis H0: pBaseline == pTrial against alternative + hypothesis H1: pBaseline != pTrial by summing p-values conditioned on individual baseline + success counts. This provides a more accurate correction for multiple testing but scales like + O(sqrt(this.baseline.numSamples)), so can eventually get slow. In that case we fall back to + zTest(). - See http://en.wikipedia.org/wiki/Statistical_hypothesis_testing#Common_test_statistics, - "Two-proportion z-test, pooled for d0 = 0". + Lower coverageAlpha increases accuracy at the cost of longer runtime. Roughly, the result + will be accurate within no more than coverageAlpha (but this ignores error due to the normal + approximation so isn't guaranteed). */ - zTest: function(zMultiplier) { - var pooledStats = new ABTest.Proportion( - this.baseline.numSuccesses + this.trial.numSuccesses, - this.baseline.numSamples + this.trial.numSamples); - var pooledPValue = pooledStats.pEstimate().value; - var pooledVarianceOfDifference = ( - pooledPValue * (1 - pooledPValue) - * (1.0 / this.baseline.numSamples + 1.0 / this.trial.numSamples)); - var pooledStandardErrorOfDifference = Math.sqrt(pooledVarianceOfDifference); - var testZValue = - Math.abs(this.differenceEstimate().value) / pooledStandardErrorOfDifference; - var adjustedOneTailedPValue = this._normal.survival(testZValue * zMultiplier); - return 2 * adjustedOneTailedPValue; + iteratedTest: function(numTrials, coverageAlpha) { + var values = this._computeTestValues(); + var trialDistribution = new ABTest.BinomialDistribution(this.trial.numSamples, + values.pooledProportion); + var baselineDistribution = new ABTest.BinomialDistribution(this.baseline.numSamples, + values.pooledProportion); + + // compute smallest and largest differences between success counts that are "at least as + // extreme" as the observed difference (the observed difference is equal to one of these) + var minExtremeDifference = + Math.floor(values.expectedDifference - values.observedAbsoluteDeviation); + var maxExtremeDifference = + Math.ceil(values.expectedDifference + values.observedAbsoluteDeviation); + + var baselineLimits = this._binomialCoverageInterval(baselineDistribution, coverageAlpha); + var pValue = 0; + for (var baselineSuccesses = baselineLimits[0]; + baselineSuccesses <= baselineLimits[1]; + baselineSuccesses++) { + // p-value of trial success counts "at least as extreme" for this particular baseline + // success count + var pValueAtBaseline = + trialDistribution.cdf(baselineSuccesses + minExtremeDifference) + + trialDistribution.survival(baselineSuccesses + maxExtremeDifference - 1); + + // this is exact because we're conditioning on the baseline count, so the multiple + // trials are independent. + var adjustedPValue = this._probabilityUnion(pValueAtBaseline, numTrials); + + var baselineProbability = baselineDistribution.mass(baselineSuccesses); + pValue += baselineProbability * adjustedPValue; + } + + // the remaining baseline values we didn't cover contribute less than coverageAlpha to the + // sum, so adding that amount gives us a conservative upper bound. + return pValue + coverageAlpha; }, }; // numTrials: number of trials to be compared to the baseline (i.e., not including the baseline) ABTest.Experiment = function(numTrials, baselineNumSuccesses, baselineNumSamples, baseAlpha) { - this._normal = new ABTest.NormalDistribution(); + normal = new ABTest.NormalDistribution(); this._baseline = new ABTest.Proportion(baselineNumSuccesses, baselineNumSamples); - var numComparisons = Math.max(1, numTrials); + this._numComparisons = Math.max(1, numTrials); // all z-values are two-tailed - var baseZCriticalValue = this._normal.inverseSurvival(baseAlpha / 2); - var alpha = baseAlpha / numComparisons // Bonferroni correction - this._zCriticalValue = this._normal.inverseSurvival(alpha / 2); + var baseZCriticalValue = normal.inverseSurvival(baseAlpha / 2); + var alpha = baseAlpha / this._numComparisons // Bonferroni correction + this._zCriticalValue = normal.inverseSurvival(alpha / 2); // to normalize for multiple testing, rather than scaling the hypothesis test's p-value, we // scale the z-value by this amount this._zMultiplier = baseZCriticalValue / this._zCriticalValue; @@ -206,7 +292,7 @@ ABTest.Experiment.prototype = { this._trialIntervalZCriticalValue), relativeImprovement: comparison.differenceRatio().valueWithInterval( this._zCriticalValue), - pValue: comparison.zTest(this._zMultiplier), + pValue: comparison.iteratedTest(this._numComparisons, 1e-5), }; }, };
159 abtest/tests/stats_test.js
 @@ -1,67 +1,136 @@ +function addToBeNearMatcher(object) { + object.addMatchers({ + // like toBeCloseTo(), but looks at absolute error rather than rounding to a fixed + // precision + toBeNear: function(expected, maxError) { + return Math.abs(this.actual - expected) < maxError; + }, + }); +} + describe('NormalDistribution', function() { - var normal = new ABTest.NormalDistribution(); + var MAX_ERROR = 1e-8; + var normal = new ABTest.NormalDistribution(1, 2); + + beforeEach(function() { + addToBeNearMatcher(this); + }); + + it('computes density', function() { + var expectedDensities = [ + [1, 0.19947114020071635], + [3, 0.12098536225957168], + [5, 0.026995483256594031], + [-1, 0.12098536225957168], + ]; + + expectedDensities.forEach(function(values) { + expect(normal.density(values[0])).toBeNear(values[1], MAX_ERROR); + }); + }); + + it('computes CDFs and survival functions', function() { + var expectedCumulativeProbabilities = [ + [1, 0.5], + [3, 0.84134474606854293], + [5, 0.97724986805182079], + [-1, 1 - 0.84134474606854293], + ]; + + expectedCumulativeProbabilities.forEach(function(values) { + expect(normal.cdf(values[0])).toBeNear(values[1], MAX_ERROR); + expect(normal.survival(values[0])).toBeNear(1 - values[1], MAX_ERROR); + }); + }); + + it('computes inverse CDFs and survival functions', function() { + var expectedValues = [ + [0.5, 1], + [0.75, 2.3489795003921632], + [0.95, 4.2897072539029448], + [0.05, 1 - 3.2897072539029448], + ]; + + expectedValues.forEach(function(values) { + expect(normal.inverseCdf(values[0])).toBeNear(values[1], MAX_ERROR); + expect(normal.inverseSurvival(values[0])).toBeNear(1 - (values[1] - 1), MAX_ERROR); + }); + }); +}); + +describe('BinomialDistribution', function() { + var binomial = new ABTest.BinomialDistribution(1000, 0.3); + var MAX_ERROR = 5e-3; beforeEach(function() { - this.addMatchers({ - // like toBeCloseTo(), but looks at absolute error rather than rounding to some precision - toBeNear: function(expected, maxError) { - return Math.abs(this.actual - expected) < maxError; - }, + addToBeNearMatcher(this); + }); + + it('computes mass', function() { + var expectedMass = [ + [300, 0.02752100382127079], + [310, 0.02152338347988187], + [340, 0.00064472915988537168], + [280, 0.01070077909763107], + ]; + + expectedMass.forEach(function(values) { + expect(binomial.mass(values[0])).toBeNear(values[1], MAX_ERROR); }); }); it('computes CDFs and survival functions', function() { - var MAX_ERROR = 7.5e-8; - - var expectedCumulativeProbabilities = {}; - expectedCumulativeProbabilities[0] = 0.5; - expectedCumulativeProbabilities[1] = 0.84134474606854293; - expectedCumulativeProbabilities[2] = 0.97724986805182079; - expectedCumulativeProbabilities[-1] = 1 - expectedCumulativeProbabilities[1]; - - for (var zValue in expectedCumulativeProbabilities) { - expect(normal.cdf(zValue)) - .toBeNear(expectedCumulativeProbabilities[zValue], MAX_ERROR); - expect(normal.survival(zValue)) - .toBeNear(1 - expectedCumulativeProbabilities[zValue], MAX_ERROR); - } + var expectedCumulativeProbabilities = [ + [300, 0.51559351981313983], + [310, 0.76630504342015282], + [340, 0.99716213728136105], + [280, 0.088579522605989086], + ]; + + expectedCumulativeProbabilities.forEach(function(values) { + expect(binomial.cdf(values[0])).toBeNear(values[1], MAX_ERROR); + expect(binomial.survival(values[0])).toBeNear(1 - values[1], MAX_ERROR); + }); }); it('computes inverse CDFs and survival functions', function() { - var MAX_ERROR = 4.5e-4; - - var expectedZValues = {}; - expectedZValues[0.5] = 0; - expectedZValues[0.75] = 0.67448975019608171; - expectedZValues[0.95] = 1.6448536269514729; - expectedZValues[0.05] = -expectedZValues[0.95]; - - for (var probability in expectedZValues) { - expect(normal.inverseCdf(probability)) - .toBeNear(expectedZValues[probability], MAX_ERROR); - expect(normal.inverseSurvival(probability)) - .toBeNear(-expectedZValues[probability], MAX_ERROR); - } + var expectedValues = [ + [0.5, 300], + [0.75, 310], + [0.95, 324], + [0.05, 276], + ]; + + expectedValues.forEach(function(values) { + expect(binomial.inverseCdf(values[0])).toBeNear(values[1], 0.5); + expect(binomial.inverseSurvival(values[0])).toBeNear(300 - (values[1] - 300), 0.5); + }); }); }); describe('Experiment', function() { - var experiment = new ABTest.Experiment(3, 20, 1000, 0.05); + var experiment = new ABTest.Experiment(3, 20, 100, 0.05); it('computes the baseline proportion', function() { var proportion = experiment.getBaselineProportion(); - expect(proportion.value).toBe(0.02); - expect(proportion.intervalWidth).toBeCloseTo(0.0074957); - expect(proportion.range().lowerBound).toBeCloseTo(0.0125043); - expect(proportion.range().upperBound).toBeCloseTo(0.0274957); + expect(proportion.value).toBe(0.2); + //expect(proportion.intervalWidth).toBeCloseTo(0.0074957); + //expect(proportion.range().lowerBound).toBeCloseTo(0.0125043); + //expect(proportion.range().upperBound).toBeCloseTo(0.0274957); }); it('computes experiment results', function() { - var results = experiment.getResults(50, 2000); - expect(results.proportion.value).toBeCloseTo(0.025); - expect(results.proportion.intervalWidth).toBeCloseTo(0.0059097); - expect(results.relativeImprovement.value).toBeCloseTo(0.25); - expect(results.relativeImprovement.intervalWidth).toBeCloseTo(0.6748677); - expect(results.pValue).toBeCloseTo(0.4838344); + var results = experiment.getResults(50, 150); + expect(results.proportion.value).toBeCloseTo(0.333333); + //expect(results.proportion.intervalWidth).toBeCloseTo(0.0059097); + expect(results.relativeImprovement.value).toBeCloseTo(0.6666667); + //expect(results.relativeImprovement.intervalWidth).toBeCloseTo(0.6748677); + expect(results.pValue).toBeCloseTo(0.0781727); + }); + + it('computes experiment results for large problems', function() { + experiment = new ABTest.Experiment(3, 50000, 100000, 0.05); + var results = experiment.getResults(101000, 200000); + expect(results.pValue).toBeCloseTo(0.0424500); }); });
1 abtest/tests/test_runner.html
 @@ -13,6 +13,7 @@ +
1 demo/abtest.html
 @@ -5,6 +5,7 @@ +
5 lib/jstat-min-08fb5e8.js
 @@ -0,0 +1,5 @@ +/** + * jStat - JavaScript Statistical Library + * Copyright (c) 2011 + * This document is licensed as free software under the terms of the + * MIT License: http://www.opensource.org/licenses/mit-license.php */this.j\$=this.jStat=function(a,b){function i(){return new i.fn.init(arguments)}var c=Array.prototype.slice,d=Object.prototype.toString,e=function(a,b){return a-b},f=function(b,c){var d=b>c?b:c;return a.pow(10,17-~~(a.log(d>0?d:-d)*a.LOG10E))},g=Array.isArray||function(a){return d.call(a)==="[object Array]"},h=function(a){return d.call(a)==="[object Function]"};i.fn=i.prototype={constructor:i,init:function(a){if(g(a[0]))if(g(a[0][0])){for(var b=0;b1?c.call(this):c.call(this)[0]},push:[].push,sort:[].sort,splice:[].splice},i.fn.init.prototype=i.fn,i.utils={calcRdx:f,isArray:g,isFunction:h},i.extend=function(a){var b=c.call(arguments),d=1,e;if(b.length===1){for(e in a)i[e]=a[e];return this}for(;d=0)d=f(b,a[c]),b=(b*d+a[c]*d)/d;return b},sumsqrd:function(a){var b=0,c=a.length;while(--c>=0)b+=a[c]*a[c];return b},sumsqerr:function(b){var c=i.mean(b),d=0,e=b.length;while(--e>=0)d+=a.pow(b[e]-c,2);return d},product:function(a){var b=1,c=a.length;while(--c>=0)b*=a[c];return b},min:function(a){var b=a[0],c=0;while(++cb&&(b=a[c]);return b},mean:function(a){return i.sum(a)/a.length},meansqerr:function(a){return i.sumsqerr(a)/a.length},geomean:function(b){return a.pow(i.product(b),1/b.length)},median:function(a){var b=a.length,c=a.slice().sort(e);return b&1?c[b/2|0]:(c[b/2-1]+c[b/2])/2},cumsum:function(a){var b=[a[0]],c=a.length,d=1;for(;df?(i=c[h],f=d,d=1,g=0):d===f?g++:d=1;return g===0?i:!1},range:function(a){var b=a.slice().sort(e);return b[b.length-1]-b[0]},variance:function(b,c){var d=i.mean(b),e=0,f=b.length-1;for(;f>=0;f--)e+=a.pow(b[f]-d,2);return e/(b.length-(c?1:0))},stdev:function(b,c){return a.sqrt(i.variance(b,c))},meandev:function(b){var c=0,d=i.mean(b),e=b.length-1;for(;e>=0;e--)c+=a.abs(b[e]-d);return c/b.length},meddev:function(b){var c=0,d=i.median(b),e=b.length-1;for(;e>=0;e--)c+=a.abs(b[e]-d);return c/b.length},quartiles:function(b){var c=b.length,d=b.slice().sort(e);return[d[a.round(c/4)-1],d[a.round(c/2)-1],d[a.round(c*3/4)-1]]},covariance:function(a,b){var c=i.mean(a),d=i.mean(b),e=[],f=a.length,g=0;for(;g1){f=b===!0?this:this.transpose();for(;e=0;a--,c++)b[c]=[this[c][a]];return i(b)},map:function(a,b){return i(i.map(this,a,b))},alter:function(a){i.alter(this,a);return this}});return i}(Math),function(a,b){(function(b){for(var c=0;c1||c<0?0:b.pow(c,d-1)*b.pow(1-c,e-1)/a.betafn(d,e)},cdf:function(b,c,d){return b>1||b<0?(b>1)*1:a.incompleteBeta(b,c,d)},inv:function(b,c,d){return a.incompleteBetaInv(b,c,d)},mean:function(a,b){return a/(a+b)},median:function(a,b){},mode:function(a,c){return a*c/(b.pow(a+c,2)*(a+c+1))},sample:function(b,c){var d=a.randg(b);return d/(d+a.randg(c))},variance:function(a,c){return a*c/(b.pow(a+c,2)*(a+c+1))}}),a.extend(a.centralF,{pdf:function(c,d,e){return c>=0?b.sqrt(b.pow(d*c,d)*b.pow(e,e)/b.pow(d*c+e,d+e))/(c*a.betafn(d/2,e/2)):undefined},cdf:function(b,c,d){return a.incompleteBeta(c*b/(c*b+d),c/2,d/2)},inv:function(b,c,d){return d/(c*(1/a.incompleteBetaInv(b,c/2,d/2)-1))},mean:function(a,b){return b>2?b/(b-2):undefined},mode:function(a,b){return a>2?b*(a-2)/(a*(b+2)):undefined},sample:function(b,c){var d=a.randg(b/2)*2,e=a.randg(c/2)*2;return d/b/(e/c)},variance:function(a,b){return b>4?2*b*b*(a+b-2)/(a*(b-2)*(b-2)*(b-4)):undefined}}),a.extend(a.cauchy,{pdf:function(a,c,d){return d/(b.pow(a-c,2)+b.pow(d,2))/b.PI},cdf:function(a,c,d){return b.atan((a-c)/d)/b.PI+.5},inv:function(a,c,d){return c+d*b.tan(b.PI*(a-.5))},median:function(a,b){return a},mode:function(a,b){return a},sample:function(c,d){return a.randn()*b.sqrt(1/(2*a.randg(.5)))*d+c}}),a.extend(a.chisquare,{pdf:function(c,d){return b.exp((d/2-1)*b.log(c)-c/2-d/2*b.log(2)-a.gammaln(d/2))},cdf:function(b,c){return a.gammap(c/2,b/2)},inv:function(b,c){return 2*a.gammapinv(b,.5*c)},mean:function(a){return a},median:function(a){return a*b.pow(1-2/(9*a),3)},mode:function(a){return a-2>0?a-2:0},sample:function(b){return a.randg(b/2)*2},variance:function(a){return 2*a}}),a.extend(a.exponential,{pdf:function(a,c){return a<0?0:c*b.exp(-c*a)},cdf:function(a,c){return a<0?0:1-b.exp(-c*a)},inv:function(a,c){return-b.log(1-a)/c},mean:function(a){return 1/a},median:function(a){return 1/a*b.log(2)},mode:function(a){return 0},sample:function(a){return-1/a*b.log(b.random())},variance:function(a){return b.pow(a,-2)}}),a.extend(a.gamma,{pdf:function(c,d,e){return b.exp((d-1)*b.log(c)-c/e-a.gammaln(d)-d*b.log(e))},cdf:function(b,c,d){return a.gammap(c,b/d)},inv:function(b,c,d){return a.gammapinv(b,c)*d},mean:function(a,b){return a*b},mode:function(a,b){if(a>1)return(a-1)*b;return undefined},sample:function(b,c){return a.randg(b)*c},variance:function(a,b){return a*b*b}}),a.extend(a.invgamma,{pdf:function(c,d,e){return b.exp(-(d+1)*b.log(c)-e/c-a.gammaln(d)+d*b.log(e))},cdf:function(b,c,d){return 1-a.gammap(c,d/b)},inv:function(b,c,d){return d/a.gammapinv(1-b,c)},mean:function(a,b){return a>1?b/(a-1):undefined},mode:function(a,b){return b/(a+1)},sample:function(b,c){return c/a.randg(b)},variance:function(a,b){return a>2?b*b/((a-1)*(a-1)*(a-2)):undefined}}),a.extend(a.kumaraswamy,{pdf:function(a,c,d){return b.exp(b.log(c)+b.log(d)+(c-1)*b.log(a)+(d-1)*b.log(1-b.pow(a,c)))},cdf:function(a,c,d){return 1-b.pow(1-b.pow(a,c),d)},mean:function(b,c){return c*a.gammafn(1+1/b)*a.gammafn(c)/a.gammafn(1+1/b+c)},median:function(a,c){return b.pow(1-b.pow(2,-1/c),1/a)},mode:function(a,c){return a>=1&&c>=1&&a!==1&&c!==1?b.pow((a-1)/(a*c-1),1/a):undefined},variance:function(a,b){}}),a.extend(a.lognormal,{pdf:function(a,c,d){return b.exp(-b.log(a)-.5*b.log(2*b.PI)-b.log(d)-b.pow(b.log(a)-c,2)/(2*d*d))},cdf:function(c,d,e){return.5+.5*a.erf((b.log(c)-d)/b.sqrt(2*e*e))},inv:function(c,d,e){return b.exp(-1.4142135623730951*e*a.erfcinv(2*c)+d)},mean:function(a,c){return b.exp(a+c*c/2)},median:function(a,c){return b.exp(a)},mode:function(a,c){return b.exp(a-c*c)},sample:function(c,d){return b.exp(a.randn()*d+c)},variance:function(a,c){return(b.exp(c*c)-1)*b.exp(2*a+c*c)}}),a.extend(a.normal,{pdf:function(a,c,d){return b.exp(-0.5*b.log(2*b.PI)-b.log(d)-b.pow(a-c,2)/(2*d*d))},cdf:function(c,d,e){return.5*(1+a.erf((c-d)/b.sqrt(2*e*e)))},inv:function(b,c,d){return-1.4142135623730951*d*a.erfcinv(2*b)+c},mean:function(a,b){return a},median:function(a,b){return a},mode:function(a,b){return a},sample:function(b,c){return a.randn()*c+b},variance:function(a,b){return b*b}}),a.extend(a.pareto,{pdf:function(a,c,d){return a>c?d*b.pow(c,d)/b.pow(a,d+1):undefined},cdf:function(a,c,d){return 1-b.pow(c/a,d)},mean:function(a,c){return c>1?c*b.pow(a,c)/(c-1):undefined},median:function(a,c){return a*c*b.SQRT2},mode:function(a,b){return a},variance:function(a,c){return c>2?a*a*c/(b.pow(c-1,2)*(c-2)):undefined}}),a.extend(a.studentt,{pdf:function(c,d){return a.gammafn((d+1)/2)/(b.sqrt(d*b.PI)*a.gammafn(d/2))*b.pow(1+c*c/d,-((d+1)/2))},cdf:function(c,d){var e=d/2;return a.incompleteBeta((c+b.sqrt(c*c+d))/(2*b.sqrt(c*c+d)),e,e)},inv:function(c,d){var e=a.incompleteBetaInv(2*b.min(c,1-c),.5*d,.5);e=b.sqrt(d*(1-e)/e);return c>0?e:-e},mean:function(a){return a>1?0:undefined},median:function(a){return 0},mode:function(a){return 0},sample:function(c){return a.randn()*b.sqrt(c/(2*a.randg(c/2)))},variance:function(a){return a>2?a/(a-2):a>1?Infinity:undefined}}),a.extend(a.weibull,{pdf:function(a,c,d){return a<0?0:d/c*b.pow(a/c,d-1)*b.exp(-b.pow(a/c,d))},cdf:function(a,c,d){return a<0?0:1-b.exp(-b.pow(a/c,d))},inv:function(a,c,d){return c*b.pow(-b.log(1-a),1/d)},mean:function(b,c){return b*a.gammafn(1+1/c)},median:function(a,c){return a*b.pow(b.log(2),1/c)},mode:function(a,c){return c>1?a*b.pow((c-1)/c,1/c):undefined},sample:function(a,c){return a*b.pow(-b.log(b.random()),1/c)},variance:function(c,d){return c*c*a.gammafn(1+2/d)-b.pow(this.mean(c,d),2)}}),a.extend(a.uniform,{pdf:function(a,b,c){return ac?0:1/(c-b)},cdf:function(a,b,c){if(ae);return d-1}}),a.extend(a.triangular,{pdf:function(a,b,c,d){return c<=b||dc?undefined:ac?0:a<=d?2*(a-b)/((c-b)*(d-b)):2*(c-a)/((c-b)*(c-d))},cdf:function(a,c,d,e){if(d<=c||ed)return undefined;if(a(a+c)/2)return a+b.sqrt((c-a)*(d-a))/b.sqrt(2)},mode:function(a,b,c){return c},sample:function(a,c,d){var e=b.random();return e<(d-a)/(c-a)?a+b.sqrt(e*(c-a)*(d-a)):c-b.sqrt((1-e)*(c-a)*(c-d))},variance:function(a,b,c){return(a*a+b*b+c*c-a*b-a*c-b*c)/18}})}(this.jStat,Math),function(a,b){function c(a,c,d){var e=1e-30,f=1,g,h,i,j,k,l,m,n,o;m=c+d,o=c+1,n=c-1,i=1,j=1-m*a/o,b.abs(j)i)for(j=0;j=1?c:1/c,p=-~(b.log(o)*8.5+c*.4+17),q,r;if(d<0||c<=0)return NaN;if(d170||d>170?b.exp(a.combinationln(c,d)):a.factorial(c)/a.factorial(d)/a.factorial(c-d)},combinationln:function(b,c){return a.factorialln(b)-a.factorialln(c)-a.factorialln(b-c)},permutation:function(b,c){return a.factorial(b)/a.factorial(b-c)},betafn:function(c,d){if(c<=0||d<=0)return undefined;return c+d>170?b.exp(a.betaln(c,d)):a.gammafn(c)*a.gammafn(d)/a.gammafn(c+d)},betaln:function(b,c){return a.gammaln(b)+a.gammaln(c)-a.gammaln(b+c)},gammapinv:function(c,d){var e=0,f=d-1,g=1e-8,h=a.gammaln(d),i,j,k,l,m,n,o;if(c>=1)return b.max(100,d+100*b.sqrt(d));if(c<=0)return 0;d>1?(n=b.log(f),o=b.exp(f*(n-1)-h),m=c<.5?c:1-c,k=b.sqrt(-2*b.log(m)),i=(2.30753+k*.27061)/(1+k*(.99229+k*.04481))-k,c<.5&&(i=-i),i=b.max(.001,d*b.pow(1-1/(9*d)-i/(3*b.sqrt(d)),3))):(k=1-d*(.253+d*.12),c1?k=o*b.exp(-(i-f)+f*(b.log(i)-n)):k=b.exp(-i+f*b.log(i)-h),l=j/k,i-=k=l/(1-.5*b.min(1,l*((d-1)/i-1))),i<=0&&(i=.5*(i+k));if(b.abs(k)0;d--)j=f,f=i*f-g+c[d],g=j;k=h*b.exp(-a*a+.5*(c[0]+i*f)-g);return e?k-1:1-k},erfc:function(b){return 1-a.erf(b)},erfcinv:function(c){var d=0,e,f,g,h;if(c>=2)return-100;if(c<=0)return 100;h=c<1?c:2-c,g=b.sqrt(-2*b.log(h/2)),e=-0.70711*((2.30753+g*.27061)/(1+g*(.99229+g*.04481))-g);for(;d<2;d++)f=a.erfc(e)-h,e+=f/(1.1283791670955126*b.exp(-e*e)-e*f);return c<1?e:-e},incompleteBetaInv:function(c,d,e){var f=1e-8,g=d-1,h=e-1,i=0,j,k,l,m,n,o,p,q,r,s,t;if(c<=0)return 0;if(c>=1)return 1;d<1||e<1?(j=b.log(d/(d+e)),k=b.log(e/(d+e)),m=b.exp(d*j)/d,n=b.exp(e*k)/e,s=m+n,c=1&&(p=.5*(p+m+1));if(b.abs(m)0)break}return p},incompleteBeta:function(d,e,f){var g=d===0||d===1?0:b.exp(a.gammaln(e+f)-a.gammaln(e)-a.gammaln(f)+e*b.log(d)+f*b.log(1-d));if(d<0||d>1)return!1;if(d<(e+1)/(e+f+2))return g*c(d,e,f)/e;return 1-g*c(1-d,f,e)/f},randn:function(c,d){var e,f,g,h,i,j;d||(d=c);if(c){j=a.zeros(c,d),j.alter(function(){return a.randn()});return j}do e=b.random(),f=1.7156*(b.random()-.5),g=e-.449871,h=b.abs(f)+.386595,i=g*g+h*(.196*h-.25472*g);while(i>.27597&&(i>.27846||f*f>-4*b.log(e)*e*e));return f/e},randg:function(c,d,e){var f=c,g,h,i,j,k,l;e||(e=d),c||(c=1);if(d){l=a.zeros(d,e),l.alter(function(){return a.randg(c)});return l}c<1&&(c+=1),g=c-1/3,h=1/b.sqrt(9*g);do{do k=a.randn(),j=1+h*k;while(j<=0);j=j*j*j,i=b.random()}while(i>1-.331*b.pow(k,4)&&b.log(i)>.5*k*k+g*(1-j+b.log(j)));if(c==f)return g*j;do i=b.random();while(i===0);return b.pow(i,1/f)*g*j}}),function(b){for(var c=0;c=0;e--){j=0;for(f=e+1;f<=g-1;f++)j=k[f]*c[e][f];k[e]=(c[e][l-1]-j)/c[e][e]}return k},gauss_jordan:function(c,d){var e=0,f=0,g=c.length,h=c[0].length,i=1,j=0,k=[],l,m,n,o;c=a.augment(c,d),n=c[0].length;for(;eh?(j[g][h]=c[g][h],k[g][h]=l[g][h]=0):gf)m=p,p=a.add(a.multiply(o,m),n),g++;return p},gauss_seidel:function(c,d,e,f){var g=0,h=c.length,i=[],j=[],k=[],l,m,n,o,p;for(;gl?(i[g][l]=c[g][l],j[g][l]=k[g][l]=0):gf)m=p,p=a.add(a.multiply(o,m),n),g=g+1;return p},SOR:function(c,d,e,f,g){var h=0,i=c.length,j=[],k=[],l=[],m,n,o,p,q;for(;hm?(j[h][m]=c[h][m],k[h][m]=l[h][m]=0):hf)n=q,q=a.add(a.multiply(p,n),o),h++;return q},householder:function(c){var d=c.length,e=c[0].length,f=0,g=[],h=[],i,j,k,l,m;for(;f0?-1:1,i=m*b.sqrt(i),j=b.sqrt((i*i-c[f+1][f]*i)/2),g=a.zeros(d,1),g[f+1][0]=(c[f+1][f]-i)/(2*j);for(k=f+2;k0?-1:1,l=o*b.sqrt(l),m=b.sqrt((l*l-c[g+1][g]*l)/2),h=a.zeros(e,1),h[g+1][0]=(c[g+1][g]-l)/(2*m);for(n=g+2;n=0;g--){p=0;for(k=g+1;k<=f-1;k++)p=j[k]*c[g][k];j[g]=d[g][0]/c[g][g]}return j},jacobi:function(c){var d=1,e=0,f=c.length,g=a.identity(f,f),h=[],i,j,k,l,m,n,o,p;while(d===1){e++,n=c[0][1],l=0,m=1;for(j=0;j0?b.PI/4:-b.PI/4:o=b.atan(2*c[l][m]/(c[l][l]-c[m][m]))/2,p=a.identity(f,f),p[l][l]=b.cos(o),p[l][m]=-b.sin(o),p[m][l]=b.sin(o),p[m][m]=b.cos(o),g=a.multiply(g,p),i=a.multiply(a.multiply(a.inv(p),c),p),c=i,d=0;for(j=1;j.001&&(d=1)}for(j=0;j=h)l=f(a,d+e),m=f(a,d),j[i]=(c[l]-2*c[m]+c[2*m-l])/(e*e),e/=2,i++;o=j.length,n=1;while(o!=1){for(p=0;pd)break;g-=1;return c[g]+(d-b[g])*m[g]+a.sq(d-b[g])*k[g]+(d-b[g])*a.sq(d-b[g])*n[g]},gauss_quadrature:function(){},PCA:function(b){var c=b.length,d=b[0].length,e=!1,f=0,g,h,i=[],j=[],k=[],l=[],m=[],n=[],o=[],p=[],q=[],r=[];for(f=0;f