--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/hexagram-6ae12361157c/hexagram/jstat-1.0.0.js	Tue Oct 22 14:17:59 2013 -0400
@@ -0,0 +1,2576 @@
+function jstat(){}
+j = jstat;
+/* Simple JavaScript Inheritance
+ * By John Resig http://ejohn.org/
+ * MIT Licensed.
+ */
+// Inspired by base2 and Prototype
+(function(){
+    var initializing = false, fnTest = /xyz/.test(function(){
+        xyz;
+    }) ? /\b_super\b/ : /.*/;
+    // The base Class implementation (does nothing)
+    this.Class = function(){};
+
+    // Create a new Class that inherits from this class
+    Class.extend = function(prop) {
+        var _super = this.prototype;
+
+        // Instantiate a base class (but only create the instance,
+        // don't run the init constructor)
+        initializing = true;
+        var prototype = new this();
+        initializing = false;
+
+        // Copy the properties over onto the new prototype
+        for (var name in prop) {
+            // Check if we're overwriting an existing function
+            prototype[name] = typeof prop[name] == "function" &&
+            typeof _super[name] == "function" && fnTest.test(prop[name]) ?
+            (function(name, fn){
+                return function() {
+                    var tmp = this._super;
+
+                    // Add a new ._super() method that is the same method
+                    // but on the super-class
+                    this._super = _super[name];
+
+                    // The method only need to be bound temporarily, so we
+                    // remove it when we're done executing
+                    var ret = fn.apply(this, arguments);
+                    this._super = tmp;
+
+                    return ret;
+                };
+            })(name, prop[name]) :
+            prop[name];
+        }
+
+        // The dummy class constructor
+        function Class() {
+            // All construction is actually done in the init method
+            if ( !initializing && this.init )
+                this.init.apply(this, arguments);
+        }
+
+        // Populate our constructed prototype object
+        Class.prototype = prototype;
+
+        // Enforce the constructor to be what we expect
+        Class.constructor = Class;
+
+        // And make this class extendable
+        Class.extend = arguments.callee;
+
+        return Class;
+    };
+})();
+
+/******************************************************************************/
+/*                           Constants                                        */
+/******************************************************************************/
+jstat.ONE_SQRT_2PI  =   0.3989422804014327;
+jstat.LN_SQRT_2PI   =   0.9189385332046727417803297;
+jstat.LN_SQRT_PId2  =   0.225791352644727432363097614947;
+jstat.DBL_MIN       =   2.22507e-308;
+jstat.DBL_EPSILON   =   2.220446049250313e-16;
+jstat.SQRT_32       =   5.656854249492380195206754896838;
+jstat.TWO_PI        =   6.283185307179586;
+jstat.DBL_MIN_EXP   =   -999;
+jstat.SQRT_2dPI     =   0.79788456080287;
+jstat.LN_SQRT_PI    =   0.5723649429247;
+/******************************************************************************/
+/*                          jstat   Functions                                 */
+/******************************************************************************/
+jstat.seq = function(min, max, length) {
+    var r = new Range(min, max, length);
+    return r.getPoints();
+}
+
+jstat.dnorm = function(x, mean, sd, log) {
+    if(mean == null) mean = 0;
+    if(sd == null) sd = 1;
+    if(log == null) log = false;
+    var n = new NormalDistribution(mean, sd);
+    if(!isNaN(x)) {
+        // is a number
+        return n._pdf(x, log);
+    } else if(x.length) {
+        var res = [];
+        for(var i = 0; i < x.length; i++) {
+            res.push(n._pdf(x[i], log));
+        }
+        return res;
+    } else {
+        throw "Illegal argument: x";
+    }
+}
+
+jstat.pnorm = function(q, mean, sd, lower_tail, log) {
+    if(mean == null) mean = 0;
+    if(sd == null) sd = 1;
+    if(lower_tail == null) lower_tail = true;
+    if(log == null) log = false;
+
+    var n = new NormalDistribution(mean, sd);
+    if(!isNaN(q)) {
+        // is a number
+        return n._cdf(q, lower_tail, log);
+    } else if(q.length) {
+        var res = [];
+        for(var i = 0; i < q.length; i++) {
+            res.push(n._cdf(q[i], lower_tail, log));
+        }
+        return res;
+    } else {
+        throw "Illegal argument: x";
+    }
+}
+
+jstat.dlnorm = function(x, meanlog, sdlog, log) {
+    if(meanlog == null) meanlog = 0;
+    if(sdlog == null) sdlog = 1;
+    if(log == null) log = false;
+    var n = new LogNormalDistribution(meanlog, sdlog);
+    if(!isNaN(x)) {
+        // is a number
+        return n._pdf(x, log);
+    } else if(x.length) {
+        var res = [];
+        for(var i = 0; i < x.length; i++) {
+            res.push(n._pdf(x[i], log));
+        }
+        return res;
+    } else {
+        throw "Illegal argument: x";
+    }
+}
+
+jstat.plnorm = function(q, meanlog, sdlog, lower_tail, log) {
+    if(meanlog == null) meanlog = 0;
+    if(sdlog == null) sdlog = 1;
+    if(lower_tail == null) lower_tail = true;
+    if(log == null) log = false;
+
+    var n = new LogNormalDistribution(meanlog, sdlog);
+    if(!isNaN(q)) {
+        // is a number
+        return n._cdf(q, lower_tail, log);
+    }
+    else if(q.length) {
+        var res = [];
+        for(var i = 0; i < q.length; i++) {
+            res.push(n._cdf(q[i], lower_tail, log));
+        }
+        return res;
+    } else {
+        throw "Illegal argument: x";
+    }
+}
+
+jstat.dbeta = function(x, alpha, beta, ncp, log) {
+    if(ncp == null) ncp = 0;
+    if(log == null) log = false;
+    var b = new BetaDistribution(alpha, beta);
+    if(!isNaN(x)) {
+        // is a number
+        return b._pdf(x, log);
+    }
+    else if(x.length) {
+        var res = [];
+        for(var i = 0; i < x.length; i++) {
+            res.push(b._pdf(x[i], log));
+        }
+        return res;
+    } else {
+        throw "Illegal argument: x";
+    }
+}
+
+jstat.pbeta = function(q, alpha, beta, ncp, lower_tail, log) {
+    if(ncp == null) ncp = 0;
+    if(log == null) log = false;
+    if(lower_tail == null) lower_tail = true;
+
+    var b = new BetaDistribution(alpha, beta);
+    if(!isNaN(q)) {
+        // is a number
+        return b._cdf(q, lower_tail, log);
+    } else if(q.length) {
+        var res = [];
+        for(var i = 0; i < q.length; i++) {
+            res.push(b._cdf(q[i], lower_tail, log));
+        }
+        return res;
+    }
+    else {
+        throw "Illegal argument: x";
+    }
+}
+
+jstat.dgamma = function(x, shape, rate, scale, log) {
+    if(rate == null) rate = 1;
+    if(scale == null) scale = 1/rate;
+    if(log == null) log = false;
+
+    var g = new GammaDistribution(shape, scale);
+    if(!isNaN(x)) {
+        // is a number
+        return g._pdf(x, log);
+    } else if(x.length) {
+        var res = [];
+        for(var i = 0; i < x.length; i++) {
+            res.push(g._pdf(x[i], log));
+        }
+        return res;
+    } else {
+        throw "Illegal argument: x";
+    }
+}
+
+jstat.pgamma = function(q, shape, rate, scale, lower_tail, log) {
+    if(rate == null) rate = 1;
+    if(scale == null) scale = 1/rate;
+    if(lower_tail == null) lower_tail = true;
+    if(log == null) log = false;
+
+    var g = new GammaDistribution(shape, scale);
+    if(!isNaN(q)) {
+        // is a number
+        return g._cdf(q, lower_tail, log);
+    } else if(q.length) {
+        var res = [];
+        for(var i = 0; i < q.length; i++) {
+            res.push(g._cdf(q[i], lower_tail, log));
+        }
+        return res;
+    } else {
+        throw "Illegal argument: x";
+    }
+
+}
+
+jstat.dt = function(x, df, ncp, log) {
+    if(log == null) log = false;
+
+    var t = new StudentTDistribution(df, ncp);
+    if(!isNaN(x)) {
+        // is a number
+        return t._pdf(x, log);
+    } else if(x.length) {
+        var res = [];
+        for(var i = 0; i < x.length; i++) {
+            res.push(t._pdf(x[i], log));
+        }
+        return res;
+    } else {
+        throw "Illegal argument: x";
+    }
+
+}
+
+jstat.pt = function(q, df, ncp, lower_tail, log) {
+    if(lower_tail == null) lower_tail = true;
+    if(log == null) log = false;
+
+    var t = new StudentTDistribution(df, ncp);
+    if(!isNaN(q)) {
+        // is a number
+        return t._cdf(q, lower_tail, log);
+    } else if(q.length) {
+        var res = [];
+        for(var i = 0; i < q.length; i++) {
+            res.push(t._cdf(q[i], lower_tail, log));
+        }
+        return res;
+    } else {
+        throw "Illegal argument: x";
+    }
+
+}
+
+jstat.plot = function(x, y, options) {
+    if(x == null) {
+        throw "x is undefined in jstat.plot";
+    }
+    if(y == null) {
+        throw "y is undefined in jstat.plot";
+    }
+    if(x.length != y.length) {
+        throw "x and y lengths differ in jstat.plot";
+    }
+
+    var flotOpt = {
+        series: {
+            lines: {
+
+            },
+            points: {
+
+        }
+        }
+    };
+
+    // combine x & y
+    var series = [];
+    if(x.length == undefined) {
+        // single point
+        series.push([x, y]);
+        flotOpt.series.points.show = true;
+    } else {
+        // array
+        for(var i = 0; i < x.length; i++) {
+            series.push([x[i], y[i]]);
+        }
+    }
+
+    var title = 'jstat graph';
+
+    // configure Flot options
+    if(options != null) {
+        // options = JSON.parse(String(options));
+        if(options.type != null) {
+            if(options.type == 'l') {
+                flotOpt.series.lines.show =  true;
+            } else if (options.type == 'p') {
+                flotOpt.series.lines.show = false;
+                flotOpt.series.points.show = true;
+            }
+        }
+        if(options.hover != null) {
+            flotOpt.grid = {
+                hoverable: options.hover
+            }
+        }
+
+        if(options.main != null) {
+            title = options.main;
+        }
+    }
+    var now = new Date();
+    var hash = now.getMilliseconds() * now.getMinutes() + now.getSeconds();
+    $('body').append('<div title="' + title + '" style="display: none;" id="'+ hash +'"><div id="graph-' + hash + '" style="width:95%; height: 95%"></div></div>');
+
+    $('#' + hash).dialog({
+        modal: false,
+        width: 475,
+        height: 475,
+        resizable: true,
+        resize: function() {
+            $.plot($('#graph-' + hash), [series], flotOpt);
+        },
+        open: function(event, ui) {
+            var id = '#graph-' + hash;
+            $.plot($('#graph-' + hash), [series], flotOpt);
+        }
+    })
+}
+
+/******************************************************************************/
+/*                          Special Functions                                 */
+/******************************************************************************/
+
+jstat.log10 = function(arg) {
+    return Math.log(arg) / Math.LN10;
+}
+
+/*
+ *
+ */
+jstat.toSigFig = function(num, n) {
+    if(num == 0) {
+        return 0;
+    }
+    var d = Math.ceil(jstat.log10(num < 0 ? -num: num));
+    var power = n - parseInt(d);
+    var magnitude = Math.pow(10,power);
+    var shifted = Math.round(num*magnitude);
+    return shifted/magnitude;
+}
+
+jstat.trunc = function(x) {
+    return (x > 0) ? Math.floor(x) : Math.ceil(x);
+}
+
+/**
+ *  Tests whether x is a finite number
+ */
+jstat.isFinite = function(x) {
+    return (!isNaN(x) && (x != Number.POSITIVE_INFINITY) && (x != Number.NEGATIVE_INFINITY));
+}
+
+/**
+ *      dopois_raw() computes the Poisson probability  lb^x exp(-lb) / x!.
+ *      This does not check that x is an integer, since dgamma() may
+ *      call this with a fractional x argument. Any necessary argument
+ *      checks should be done in the calling function.
+ */
+jstat.dopois_raw = function(x, lambda, give_log) {
+    /*       x >= 0 ; integer for dpois(), but not e.g. for pgamma()!
+        lambda >= 0
+     */
+    if (lambda == 0) {
+        if(x == 0) {
+            return(give_log) ? 0.0 : 1.0; //R_D__1
+        }
+        return (give_log) ? Number.NEGATIVE_INFINITY : 0.0; // R_D__0
+    }
+    if (!jstat.isFinite(lambda)) return (give_log) ? Number.NEGATIVE_INFINITY : 0.0; //R_D__0;
+    if (x < 0) return(give_log) ? Number.NEGATIVE_INFINITY : 0.0; //R_D__0
+    if (x <= lambda * jstat.DBL_MIN) {
+        return (give_log) ? -lambda : Math.exp(-lambda);    // R_D_exp(-lambda)
+    }
+    if (lambda < x * jstat.DBL_MIN) {
+        var param = -lambda + x*Math.log(lambda) -jstat.lgamma(x+1);
+        return (give_log) ? param : Math.exp(param);    // R_D_exp(-lambda + x*log(lambda) -lgammafn(x+1))
+    }
+    var param1 = jstat.TWO_PI * x;  // f
+    var param2 = -jstat.stirlerr(x)-jstat.bd0(x,lambda);    // x
+    return (give_log) ? -0.5*Math.log(param1)+param2 : Math.exp(param2)/Math.sqrt(param1);  // R_D_fexp(M_2PI*x, -stirlerr(x)-bd0(x,lambda))
+//return(R_D_fexp( , -stirlerr(x)-bd0(x,lambda) ));
+}
+
+/**	Evaluates the "deviance part"
+ *	bd0(x,M) :=  M * D0(x/M) = M*[ x/M * log(x/M) + 1 - (x/M) ] =
+ *		  =  x * log(x/M) + M - x
+ *	where M = E[X] = n*p (or = lambda), for	  x, M > 0
+ *
+ *	in a manner that should be stable (with small relative error)
+ *	for all x and M=np. In particular for x/np close to 1, direct
+ *	evaluation fails, and evaluation is based on the Taylor series
+ *	of log((1+v)/(1-v)) with v = (x-np)/(x+np).
+ */
+jstat.bd0 = function(x, np) {
+    var ej, s, s1, v, j;
+    if(!jstat.isFinite(x) || !jstat.isFinite(np) || np == 0.0) throw "illegal parameter in jstat.bd0";
+
+    if(Math.abs(x-np) > 0.1*(x+np)) {
+        v = (x-np)/(x+np);
+        s = (x-np)*v;/* s using v -- change by MM */
+        ej = 2*x*v;
+        v = v*v;
+        for (j=1; ; j++) { /* Taylor series */
+            ej *= v;
+            s1 = s+ej/((j<<1)+1);
+            if (s1==s) /* last term was effectively 0 */
+                return(s1);
+            s = s1;
+        }
+    }
+    /* else:  | x - np |  is not too small */
+    return(x*Math.log(x/np)+np-x);
+}
+
+/**    Computes the log of the error term in Stirling's formula.
+ *      For n > 15, uses the series 1/12n - 1/360n^3 + ...
+ *      For n <=15, integers or half-integers, uses stored values.
+ *      For other n < 15, uses lgamma directly (don't use this to
+ *        write lgamma!)
+ */
+jstat.stirlerr= function(n) {
+    var S0 = 0.083333333333333333333;
+    var S1 = 0.00277777777777777777778;
+    var S2 = 0.00079365079365079365079365;
+    var S3 = 0.000595238095238095238095238;
+    var S4 = 0.0008417508417508417508417508;
+
+    var sferr_halves = [
+    0.0, /* n=0 - wrong, place holder only */
+    0.1534264097200273452913848,  /* 0.5 */
+    0.0810614667953272582196702,  /* 1.0 */
+    0.0548141210519176538961390,  /* 1.5 */
+    0.0413406959554092940938221,  /* 2.0 */
+    0.03316287351993628748511048, /* 2.5 */
+    0.02767792568499833914878929, /* 3.0 */
+    0.02374616365629749597132920, /* 3.5 */
+    0.02079067210376509311152277, /* 4.0 */
+    0.01848845053267318523077934, /* 4.5 */
+    0.01664469118982119216319487, /* 5.0 */
+    0.01513497322191737887351255, /* 5.5 */
+    0.01387612882307074799874573, /* 6.0 */
+    0.01281046524292022692424986, /* 6.5 */
+    0.01189670994589177009505572, /* 7.0 */
+    0.01110455975820691732662991, /* 7.5 */
+    0.010411265261972096497478567, /* 8.0 */
+    0.009799416126158803298389475, /* 8.5 */
+    0.009255462182712732917728637, /* 9.0 */
+    0.008768700134139385462952823, /* 9.5 */
+    0.008330563433362871256469318, /* 10.0 */
+    0.007934114564314020547248100, /* 10.5 */
+    0.007573675487951840794972024, /* 11.0 */
+    0.007244554301320383179543912, /* 11.5 */
+    0.006942840107209529865664152, /* 12.0 */
+    0.006665247032707682442354394, /* 12.5 */
+    0.006408994188004207068439631, /* 13.0 */
+    0.006171712263039457647532867, /* 13.5 */
+    0.005951370112758847735624416, /* 14.0 */
+    0.005746216513010115682023589, /* 14.5 */
+    0.005554733551962801371038690  /* 15.0 */
+    ];
+
+    var nn;
+
+    if (n <= 15.0) {
+        nn = n + n;
+        if (nn == parseInt(nn)) return(sferr_halves[parseInt(nn)]);
+        return(jstat.lgamma(n + 1.0) - (n + 0.5)*Math.log(n) + n - jstat.LN_SQRT_2PI);
+    }
+
+    nn = n*n;
+    if (n>500) return((S0-S1/nn)/n);
+    if (n> 80) return((S0-(S1-S2/nn)/nn)/n);
+    if (n> 35) return((S0-(S1-(S2-S3/nn)/nn)/nn)/n);
+    /* 15 < n <= 35 : */
+    return((S0-(S1-(S2-(S3-S4/nn)/nn)/nn)/nn)/n);
+}
+
+
+
+/**    The function lgamma computes log|gamma(x)|.  The function
+ *    lgammafn_sign in addition assigns the sign of the gamma function
+ *    to the address in the second argument if this is not null.
+ */
+jstat.lgamma = function(x) {
+    function lgammafn_sign(x, sgn) {
+        var ans, y, sinpiy;
+        var xmax = 2.5327372760800758e+305;
+        var dxrel = 1.490116119384765696e-8;
+
+        // if (xmax == 0) {/* initialize machine dependent constants _ONCE_ */
+        //     xmax = jstat.DBL_MAX/Math.log(jstat.DBL_MAX);/* = 2.533 e305	 for IEEE double */
+        //     dxrel = Math.sqrt(jstat.DBL_EPSILON);/* sqrt(Eps) ~ 1.49 e-8  for IEEE double */
+        // }
+
+        /* For IEEE double precision DBL_EPSILON = 2^-52 = 2.220446049250313e-16 :
+           xmax  = DBL_MAX / log(DBL_MAX) = 2^1024 / (1024 * log(2)) = 2^1014 / log(2)
+           dxrel = sqrt(DBL_EPSILON) = 2^-26 = 5^26 * 1e-26 (is *exact* below !)
+         */
+
+        if (sgn != null) sgn = 1;
+
+        if(isNaN(x)) return x;
+
+        if (x < 0 && (Math.floor(-x) % 2.0) == 0)
+            if (sgn != null) sgn = -1;
+
+        if (x <= 0 && x == jstat.trunc(x)) { /* Negative integer argument */
+            console.warn("Negative integer argument in lgammafn_sign");
+            return Number.POSITIVE_INFINITY;/* +Inf, since lgamma(x) = log|gamma(x)| */
+        }
+
+        y = Math.abs(x);
+
+        if(y <= 10) return Math.log(Math.abs(jstat.gamma(x)));  // TODO: implement jstat.gamma
+
+        if(y > xmax) {
+            console.warn("Illegal arguement passed to lgammafn_sign");
+            return Number.POSITIVE_INFINITY;
+        }
+
+        if(x > 0) {
+            if(x > 1e17) {
+                return (x*(Math.log(x)-1.0));
+            } else if(x > 4934720.0) {
+                return (jstat.LN_SQRT_2PI + (x-0.5) * Math.log(x) - x);
+            } else {
+                return jstat.LN_SQRT_2PI + (x-0.5) * Math.log(x) - x + jstat.lgammacor(x);  // TODO: implement lgammacor
+            }
+        }
+
+        sinpiy = Math.abs(Math.sin(Math.PI * y));
+
+        if(sinpiy == 0) {
+            throw "Should never happen!!";
+        }
+
+        ans = jstat.LN_SQRT_PId2 + (x - 0.5) * Math.log(y) - x - Math.log(sinpiy) - jstat.lgammacor(y);
+
+        if(Math.abs((x-jstat.trunc(x-0.5))* ans / x) < dxrel) {
+            throw "The answer is less than half the precision argument too close to a negative integer";
+        }
+        return ans;
+    }
+
+    return lgammafn_sign(x, null);
+}
+
+jstat.gamma = function(x) {
+    var xbig = 171.624;
+    var p = [
+    -1.71618513886549492533811,
+    24.7656508055759199108314,-379.804256470945635097577,
+    629.331155312818442661052,866.966202790413211295064,
+    -31451.2729688483675254357,-36144.4134186911729807069,
+    66456.1438202405440627855
+    ];
+    var q = [
+    -30.8402300119738975254353,
+    315.350626979604161529144,-1015.15636749021914166146,
+    -3107.77167157231109440444,22538.1184209801510330112,
+    4755.84627752788110767815,-134659.959864969306392456,
+    -115132.259675553483497211
+    ];
+    var c = [
+    -.001910444077728,8.4171387781295e-4,
+    -5.952379913043012e-4,7.93650793500350248e-4,
+    -.002777777777777681622553,.08333333333333333331554247,
+    .0057083835261
+    ];
+
+    var i,n,parity,fact,xden,xnum,y,z,yi,res,sum,ysq;
+
+    parity = (0);
+    fact = 1.0;
+    n = 0;
+    y=x;
+    if(y <= 0.0) {
+        /* -------------------------------------------------------------
+	   Argument is negative
+	   ------------------------------------------------------------- */
+        y = -x;
+        yi = jstat.trunc(y);
+        res = y - yi;
+        if (res != 0.0) {
+            if (yi != jstat.trunc(yi * 0.5) * 2.0)
+                parity = (1);
+            fact = -Math.PI / Math.sin(Math.PI * res);
+            y += 1.0;
+        } else {
+            return(Number.POSITIVE_INFINITY);
+        }
+    }
+    /* -----------------------------------------------------------------
+       Argument is positive
+       -----------------------------------------------------------------*/
+    if (y < jstat.DBL_EPSILON) {
+        /* --------------------------------------------------------------
+	   Argument < EPS
+	   -------------------------------------------------------------- */
+        if (y >= jstat.DBL_MIN) {
+            res = 1.0 / y;
+        } else {
+            return(Number.POSITIVE_INFINITY);
+        }
+    } else if (y < 12.0) {
+        yi = y;
+        if (y < 1.0) {
+            /* ---------------------------------------------------------
+	       EPS < argument < 1
+	       --------------------------------------------------------- */
+            z = y;
+            y += 1.0;
+        } else {
+            /* -----------------------------------------------------------
+	       1 <= argument < 12, reduce argument if necessary
+	       ----------------------------------------------------------- */
+            n = parseInt(y) - 1;
+            y -= parseFloat(n);
+            z = y - 1.0;
+        }
+        /* ---------------------------------------------------------
+	   Evaluate approximation for 1. < argument < 2.
+	   ---------------------------------------------------------*/
+        xnum = 0.0;
+        xden = 1.0;
+        for (i = 0; i < 8; ++i) {
+            xnum = (xnum + p[i]) * z;
+            xden = xden * z + q[i];
+        }
+        res = xnum / xden + 1.0;
+        if (yi < y) {
+            /* --------------------------------------------------------
+	       Adjust result for case  0. < argument < 1.
+	       -------------------------------------------------------- */
+            res /= yi;
+        } else if (yi > y) {
+            /* ----------------------------------------------------------
+	       Adjust result for case  2. < argument < 12.
+	       ---------------------------------------------------------- */
+            for (i = 0; i < n; ++i) {
+                res *= y;
+                y += 1.0;
+            }
+        }
+    } else {
+        /* -------------------------------------------------------------
+	   Evaluate for argument >= 12.,
+	   ------------------------------------------------------------- */
+        if (y <= xbig) {
+            ysq = y * y;
+            sum = c[6];
+            for (i = 0; i < 6; ++i) {
+                sum = sum / ysq + c[i];
+            }
+            sum = sum / y - y + jstat.LN_SQRT_2PI;
+            sum += (y - 0.5) * Math.log(y);
+            res = Math.exp(sum);
+        } else {
+            return(Number.POSITIVE_INFINITY);
+        }
+    }
+    /* ----------------------------------------------------------------------
+       Final adjustments and return
+       ----------------------------------------------------------------------*/
+    if (parity)
+        res = -res;
+    if (fact != 1.0)
+        res = fact / res;
+    return res;
+}
+
+/**    Compute the log gamma correction factor for x >= 10 so that
+ *
+ *    log(gamma(x)) = .5*log(2*pi) + (x-.5)*log(x) -x + lgammacor(x)
+ *
+ *    [ lgammacor(x) is called	Del(x)	in other contexts (e.g. dcdflib)]
+ */
+jstat.lgammacor = function(x) {
+    var algmcs = [
+    +.1666389480451863247205729650822e+0,
+    -.1384948176067563840732986059135e-4,
+    +.9810825646924729426157171547487e-8,
+    -.1809129475572494194263306266719e-10,
+    +.6221098041892605227126015543416e-13,
+    -.3399615005417721944303330599666e-15,
+    +.2683181998482698748957538846666e-17,
+    -.2868042435334643284144622399999e-19,
+    +.3962837061046434803679306666666e-21,
+    -.6831888753985766870111999999999e-23,
+    +.1429227355942498147573333333333e-24,
+    -.3547598158101070547199999999999e-26,
+    +.1025680058010470912000000000000e-27,
+    -.3401102254316748799999999999999e-29,
+    +.1276642195630062933333333333333e-30
+    ];
+
+    var tmp;
+    var nalgm = 5;
+    var xbig = 94906265.62425156;
+    var xmax = 3.745194030963158e306;
+
+    if(x < 10) {
+        return Number.NaN;
+    } else if (x >= xmax) {
+        throw "Underflow error in lgammacor";
+    } else if (x < xbig) {
+        tmp = 10 / x;
+        return jstat.chebyshev(tmp*tmp*2-1,algmcs,nalgm) / x;
+    }
+    return 1 / (x*12);
+}
+
+/*
+ * Incomplete Beta function
+ */
+jstat.incompleteBeta = function(a, b, x) {
+    /*
+     * Used by incompleteBeta: Evaluates continued fraction for incomplete
+     * beta function by modified Lentz's method.
+     */
+    function betacf(a, b, x) {
+        var MAXIT = 100;
+        var EPS = 3.0e-12;
+        var FPMIN = 1.0e-30;
+
+        var m,m2,aa,c,d,del,h,qab,qam,qap;
+
+        qab=a+b;
+        qap=a+1.0;
+        qam=a-1.0;
+        c=1.0;
+        d=1.0-qab*x/qap;
+
+        if(Math.abs(d) < FPMIN) {
+            d=FPMIN;
+        }
+
+        d = 1.0/d;
+        h=d;
+        for(m = 1; m <= MAXIT; m++) {
+            m2=2*m;
+            aa=m*(b-m)*x/((qam+m2)*(a+m2));
+            d=1.0+aa*d;
+            if(Math.abs(d) < FPMIN) {
+                d = FPMIN;
+            }
+            c=1.0+aa/c;
+            if(Math.abs(c) < FPMIN) {
+                c = FPMIN;
+            }
+            d=1.0/d;
+            h *= d*c;
+            aa = -(a+m)*(qab+m)*x/((a+m2) * (qap+m2));
+            d=1.0+aa*d;
+            if(Math.abs(d) < FPMIN) {
+                d = FPMIN;
+            }
+            c=1.0+aa/c;
+            if(Math.abs(c) < FPMIN) {
+                c=FPMIN;
+            }
+            d=1.0/d;
+            del=d*c;
+            h *= del;
+            if(Math.abs(del-1.0) < EPS) {
+                // are we done?
+                break;
+            }
+        }
+        if(m > MAXIT) {
+            console.warn("a or b too big, or MAXIT too small in betacf: " + a + ", " + b + ", " + x + ", " + h);
+            return h;
+        }
+        if(isNaN(h)) {
+            console.warn(a + ", " + b + ", " + x);
+        }
+        return h;
+    }
+
+    var bt;
+
+    if(x < 0.0 || x > 1.0) {
+        throw "bad x in routine incompleteBeta";
+    }
+    if(x == 0.0 || x == 1.0) {
+        bt = 0.0;
+    } else {
+        bt = Math.exp(jstat.lgamma(a+b) - jstat.lgamma(a) - jstat.lgamma(b) + a * Math.log(x)+ b * Math.log(1.0-x));
+    }
+    if(x < (a + 1.0)/(a+b+2.0)) {
+        return bt * betacf(a,b,x)/a;
+    } else {
+        return 1.0-bt*betacf(b,a,1.0-x)/b;
+    }
+}
+
+/**   Evaluates the n-term Chebyshev series
+ *    "a" at "x".
+ */
+jstat.chebyshev = function(x, a, n) {
+    var b0, b1, b2, twox;
+    var i;
+
+    if (n < 1 || n > 1000) return Number.NaN;
+
+    if (x < -1.1 || x > 1.1) return Number.NaN;
+
+    twox = x * 2;
+    b2 = b1 = 0;
+    b0 = 0;
+    for (i = 1; i <= n; i++) {
+        b2 = b1;
+        b1 = b0;
+        b0 = twox * b1 - b2 + a[n - i];
+    }
+    return (b0 - b2) * 0.5;
+}
+
+jstat.fmin2 = function(x, y) {
+    return (x < y) ? x : y;
+}
+
+jstat.log1p = function(x) {
+    // http://kevin.vanzonneveld.net
+    // +   original by: Brett Zamir (http://brett-zamir.me)
+    // %          note 1: Precision 'n' can be adjusted as desired
+    // *     example 1: log1p(1e-15);
+    // *     returns 1: 9.999999999999995e-16
+
+    var ret = 0,
+    n = 50; // degree of precision
+    if (x <= -1) {
+        return Number.NEGATIVE_INFINITY; // JavaScript style would be to return Number.NEGATIVE_INFINITY
+    }
+    if (x < 0 || x > 1) {
+        return Math.log(1 + x);
+    }
+    for (var i = 1; i < n; i++) {
+        if ((i % 2) === 0) {
+            ret -= Math.pow(x, i) / i;
+        } else {
+            ret += Math.pow(x, i) / i;
+        }
+    }
+    return ret;
+}
+
+jstat.expm1 = function(x) {
+    var y, a  = Math.abs(x);
+    if(a < jstat.DBL_EPSILON) return x;
+    if(a > 0.697) return Math.exp(x) - 1; /* negligable cancellation */
+    if(a > 1e-8) {
+        y = Math.exp(x) - 1;
+    } else {
+        y = (x / 2 + 1) * x;
+    }
+
+    /* Newton step for solving   log(1 + y) = x   for y : */
+    /* WARNING: does not work for y ~ -1: bug in 1.5.0 */
+    y -= (1 + y) * (jstat.log1p(y) - x);
+    return y;
+}
+
+jstat.logBeta = function(a, b) {
+    var corr, p, q;
+    p = q = a;
+    if(b < p) p = b;/* := min(a,b) */
+    if(b > q) q = b;/* := max(a,b) */
+
+    /* both arguments must be >= 0 */
+    if (p < 0) {
+        console.warn('Both arguements must be >= 0');
+        return Number.NaN;
+    }
+    else if (p == 0) {
+        return Number.POSITIVE_INFINITY;
+    }
+    else if (!jstat.isFinite(q)) { /* q == +Inf */
+        return Number.NEGATIVE_INFINITY;
+    }
+
+    if (p >= 10) {
+        /* p and q are big. */
+        corr = jstat.lgammacor(p) + jstat.lgammacor(q) - jstat.lgammacor(p + q);
+        return Math.log(q) * -0.5 + jstat.LN_SQRT_2PI + corr
+        + (p - 0.5) * Math.log(p / (p + q)) + q * jstat.log1p(-p / (p + q));
+    }
+    else if (q >= 10) {
+        /* p is small, but q is big. */
+        corr = jstat.lgammacor(q) - jstat.lgammacor(p + q);
+        return jstat.lgamma(p) + corr + p - p * Math.log(p + q)
+        + (q - 0.5) * jstat.log1p(-p / (p + q));
+    }
+    else
+        /* p and q are small: p <= q < 10. */
+        return Math.log(jstat.gamma(p) * (jstat.gamma(q) / jstat.gamma(p + q)));
+}
+
+jstat.dbinom_raw = function(x, n, p, q, give_log) {
+    if(give_log == null) give_log = false;
+    var lf, lc;
+
+    if(p == 0) {
+        if(x == 0) {
+            // R_D__1
+            return (give_log) ? 0.0 : 1.0;
+        } else {
+            // R_D__0
+            return (give_log) ? Number.NEGATIVE_INFINITY : 0.0;
+        }
+    }
+    if(q == 0) {
+        if(x == n) {
+            // R_D__1
+            return (give_log) ? 0.0 : 1.0;
+        } else {
+            // R_D__0
+            return (give_log) ? Number.NEGATIVE_INFINITY : 0.0;
+        }
+    }
+
+    if (x == 0) {
+        if(n == 0) return (give_log) ? 0.0 : 1.0;   //R_D__1;
+        lc = (p < 0.1) ? -jstat.bd0(n,n*q) - n*p : n*Math.log(q);
+        return ( give_log ) ? lc : Math.exp(lc); //R_D_exp(lc)
+    }
+
+    if (x == n) {
+        lc = (q < 0.1) ? -jstat.bd0(n,n*p) - n*q : n*Math.log(p);
+        return ( give_log ) ? lc : Math.exp(lc); //R_D_exp(lc)
+    }
+
+    if (x < 0 || x > n) return (give_log) ? Number.NEGATIVE_INFINITY : 0.0; // R_D__0;
+
+    /* n*p or n*q can underflow to zero if n and p or q are small.  This
+       used to occur in dbeta, and gives NaN as from R 2.3.0.  */
+    lc = jstat.stirlerr(n) - jstat.stirlerr(x) - jstat.stirlerr(n-x) - jstat.bd0(x,n*p) - jstat.bd0(n-x,n*q);
+
+    /* f = (M_2PI*x*(n-x))/n; could overflow or underflow */
+    /* Upto R 2.7.1:
+     * lf = log(M_2PI) + log(x) + log(n-x) - log(n);
+     * -- following is much better for  x << n : */
+    lf = Math.log(jstat.TWO_PI) + Math.log(x) + jstat.log1p(- x/n);
+
+    return (give_log) ? lc - 0.5*lf : Math.exp(lc - 0.5*lf); // R_D_exp(lc - 0.5*lf);
+}
+
+jstat.max = function(values) {
+    var max = Number.NEGATIVE_INFINITY;
+    for(var i = 0; i < values.length; i++) {
+        if(values[i] > max) {
+            max = values[i];
+        }
+    }
+    return max;
+}
+
+/******************************************************************************/
+/*                      Probability Distributions                             */
+/******************************************************************************/
+
+/**
+ * Range class
+ */
+var Range = Class.extend({
+    init: function(min, max, numPoints) {
+        this._minimum = parseFloat(min);
+        this._maximum = parseFloat(max);
+        this._numPoints = parseFloat(numPoints);
+    },
+    getMinimum: function() {
+        return this._minimum;
+    },
+    getMaximum: function() {
+        return this._maximum;
+    },
+    getNumPoints: function() {
+        return this._numPoints;
+    },
+    getPoints: function() {
+        var results = [];
+        var x = this._minimum;
+        var step = (this._maximum-this._minimum)/(this._numPoints-1);
+        for(var i = 0; i < this._numPoints; i++) {
+            results[i] = parseFloat(x.toFixed(6));
+            x += step;
+        }
+        return results;
+    }
+});
+
+Range.validate = function(range) {
+    if( ! range instanceof Range) {
+        return false;
+    }
+    if(isNaN(range.getMinimum()) || isNaN(range.getMaximum()) || isNaN(range.getNumPoints()) || range.getMaximum() < range.getMinimum() || range.getNumPoints() <= 0) {
+        return false;
+    }
+    return true;
+}
+
+var ContinuousDistribution = Class.extend({
+    init: function(name) {
+        this._name = name;
+    },
+    toString: function() {
+        return this._string;
+    },
+    getName: function() {
+        return this._name;
+    },
+    getClassName: function() {
+        return this._name + 'Distribution';
+    },
+    density: function(valueOrRange) {
+        if(!isNaN(valueOrRange)) {
+            // single value
+            return parseFloat(this._pdf(valueOrRange).toFixed(15));
+        } else if (Range.validate(valueOrRange)) {
+            // multiple values
+            var points = valueOrRange.getPoints();
+
+            var result = [];
+            // For each point in the range
+            for(var i = 0; i < points.length; i++) {
+                result[i] = parseFloat(this._pdf(points[i]));
+            }
+            return result;
+        } else {
+            // neither value or range
+            throw "Invalid parameter supplied to " + this.getClassName() + ".density()";
+        }
+    },
+    cumulativeDensity: function(valueOrRange) {
+        if(!isNaN(valueOrRange)) {
+            // single value
+            return parseFloat(this._cdf(valueOrRange).toFixed(15));
+        } else if (Range.validate(valueOrRange)) {
+            // multiple values
+            var points = valueOrRange.getPoints();
+            var result = [];
+            // For each point in the range
+            for(var i = 0; i < points.length; i++) {
+                result[i] = parseFloat(this._cdf(points[i]));
+            }
+            return result;
+        } else {
+            // neither value or range
+            throw "Invalid parameter supplied to " + this.getClassName() + ".cumulativeDensity()";
+        }
+    },
+    getRange: function(standardDeviations, numPoints) {
+        if(standardDeviations == null) {
+            standardDeviations = 5;
+        }
+        if(numPoints == null) {
+            numPoints = 100;
+        }
+        var min = this.getMean() - standardDeviations * Math.sqrt(this.getVariance());
+        var max = this.getMean() + standardDeviations * Math.sqrt(this.getVariance());
+
+        if(this.getClassName() == 'GammaDistribution' || this.getClassName() == 'LogNormalDistribution') {
+            min = 0.0;
+            max = this.getMean() + standardDeviations * Math.sqrt(this.getVariance());
+        } else if(this.getClassName() == 'BetaDistribution') {
+            min = 0.0;
+            max = 1.0;
+        }
+
+
+        var range = new Range(min, max, numPoints);
+        return range;
+    },
+    getVariance: function(){},
+    getMean: function(){},
+    getQuantile: function(p) {
+        var self = this;
+        /*
+         *  Recursive function to find the closest match
+         */
+        function findClosestMatch(range, p) {
+            var ERR = 1.0e-5;
+            var xs = range.getPoints();
+            var closestIndex = 0;
+            var closestDistance = 999;
+
+            for(var i=0; i<xs.length; i++) {
+                var pp = self.cumulativeDensity(xs[i]);
+                var distance = Math.abs(pp - p);
+                if(distance < closestDistance) {
+                    // closer value found
+                    closestIndex = i;
+                    closestDistance = distance;
+                }
+            }
+            if(closestDistance <= ERR) {
+                // Acceptable - return value;
+                return xs[closestIndex];
+            } else {
+                // Calculate the new range
+                var newRange = new Range(xs[closestIndex-1], xs[closestIndex+1],20);
+                return findClosestMatch(newRange, p);
+            }
+        }
+        var range = this.getRange(5, 20);
+        return findClosestMatch(range, p);
+    }
+});
+
+/**
+ * A normal distribution object
+ */
+var NormalDistribution = ContinuousDistribution.extend({
+    init: function(mean, sigma) {
+        this._super('Normal');
+        this._mean = parseFloat(mean);
+        this._sigma = parseFloat(sigma);
+        this._string = "Normal ("+this._mean.toFixed(2)+", " + this._sigma.toFixed(2) + ")";
+    },
+    _pdf: function(x, give_log) {
+        if(give_log == null) {
+            give_log=false;
+        }  // default is false;
+        var sigma = this._sigma;
+        var mu = this._mean;
+        if(!jstat.isFinite(sigma)) {
+            return (give_log) ? Number.NEGATIVE_INFINITY : 0.0
+        }
+        if(!jstat.isFinite(x) && mu == x) {
+            return Number.NaN;
+        }
+        if(sigma<=0) {
+            if(sigma < 0) {
+                throw "invalid sigma in _pdf";
+            }
+            return (x==mu)?Number.POSITIVE_INFINITY:(give_log)?Number.NEGATIVE_INFINITY:0.0;
+        }
+        x=(x-mu)/sigma;
+        if(!jstat.isFinite(x)){
+            return (give_log)?Number.NEGATIVE_INFINITY:0.0;
+        }
+        return (give_log ? -(jstat.LN_SQRT_2PI + 0.5 * x * x + Math.log(sigma)) :
+            jstat.ONE_SQRT_2PI * Math.exp(-0.5 * x * x) / sigma);
+    },
+    _cdf: function(x, lower_tail, log_p) {
+
+        if(lower_tail == null) lower_tail = true;
+        if(log_p == null) log_p = false;
+
+        function pnorm_both(x, cum, ccum, i_tail, log_p) {
+            /*  i_tail in {0,1,2} means: "lower", "upper", or "both" :
+                if(lower) return  *cum := P[X <= x]
+                if(upper) return *ccum := P[X >  x] = 1 - P[X <= x]
+             */
+
+            var a = [
+            2.2352520354606839287,
+            161.02823106855587881,
+            1067.6894854603709582,
+            18154.981253343561249,
+            0.065682337918207449113
+            ];
+            var b = [
+            47.20258190468824187,
+            976.09855173777669322,
+            10260.932208618978205,
+            45507.789335026729956
+            ];
+            var c = [
+            0.39894151208813466764,
+            8.8831497943883759412,
+            93.506656132177855979,
+            597.27027639480026226,
+            2494.5375852903726711,
+            6848.1904505362823326,
+            11602.651437647350124,
+            9842.7148383839780218,
+            1.0765576773720192317e-8
+            ];
+            var d = [
+            22.266688044328115691,
+            235.38790178262499861,
+            1519.377599407554805,
+            6485.558298266760755,
+            18615.571640885098091,
+            34900.952721145977266,
+            38912.003286093271411,
+            19685.429676859990727
+            ];
+            var p = [
+            0.21589853405795699,
+            0.1274011611602473639,
+            0.022235277870649807,
+            0.001421619193227893466,
+            2.9112874951168792e-5,
+            0.02307344176494017303
+            ];
+            var q = [
+            1.28426009614491121,
+            0.468238212480865118,
+            0.0659881378689285515,
+            0.00378239633202758244,
+            7.29751555083966205e-5
+            ];
+
+            var xden, xnum, temp, del, eps, xsq, y, i, lower, upper;
+
+            /* Consider changing these : */
+            eps = jstat.DBL_EPSILON * 0.5;
+
+            /* i_tail in {0,1,2} =^= {lower, upper, both} */
+            lower = i_tail != 1;
+            upper = i_tail != 0;
+
+            y = Math.abs(x);
+
+            if (y <= 0.67448975) { /* qnorm(3/4) = .6744.... -- earlier had 0.66291 */
+                if (y > eps) {
+                    xsq = x * x;
+                    xnum = a[4] * xsq;
+                    xden = xsq;
+                    for (i = 0; i < 3; ++i) {
+                        xnum = (xnum + a[i]) * xsq;
+                        xden = (xden + b[i]) * xsq;
+                    }
+                } else {
+                    xnum = xden = 0.0;
+                }
+                temp = x * (xnum + a[3]) / (xden + b[3]);
+                if(lower)  cum = 0.5 + temp;
+                if(upper) ccum = 0.5 - temp;
+                if(log_p) {
+                    if(lower)  cum = Math.log(cum);
+                    if(upper) ccum = Math.log(ccum);
+                }
+
+            } else if (y <= jstat.SQRT_32) {
+                /* Evaluate pnorm for 0.674.. = qnorm(3/4) < |x| <= sqrt(32) ~= 5.657 */
+
+                xnum = c[8] * y;
+                xden = y;
+                for (i = 0; i < 7; ++i) {
+                    xnum = (xnum + c[i]) * y;
+                    xden = (xden + d[i]) * y;
+                }
+                temp = (xnum + c[7]) / (xden + d[7]);
+
+                /* do_del */
+                xsq = jstat.trunc(x * 16) / 16;
+                del = (x - xsq) * (x + xsq);
+                if(log_p) {
+                    cum = (-xsq * xsq * 0.5) + (-del * 0.5) + Math.log(temp);
+                    if((lower && x > 0.) || (upper && x <= 0.))
+                        ccum = jstat.log1p(-Math.exp(-xsq * xsq * 0.5) *
+                            Math.exp(-del * 0.5) * temp);
+                }
+                else {
+                    cum = Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp;
+                    ccum = 1.0 - cum;
+                }
+                /* end do_del */
+
+                /* swap_tail */
+                if (x > 0.0) {/* swap  ccum <--> cum */
+                    temp = cum;
+                    if(lower) {
+                        cum = ccum;
+
+                    }
+                    ccum = temp;
+                }
+            /* end swap_tail */
+
+            }
+            /* else	  |x| > sqrt(32) = 5.657 :
+             * the next two case differentiations were really for lower=T, log=F
+             * Particularly	 *not*	for  log_p !
+
+             * Cody had (-37.5193 < x  &&  x < 8.2924) ; R originally had y < 50
+             *
+             * Note that we do want symmetry(0), lower/upper -> hence use y
+             */
+
+            else if((log_p && y < 1e170)|| (lower && -37.5193 < x  &&  x < 8.2924)
+                || (upper && -8.2924  < x  &&  x < 37.5193)) {
+                /* Evaluate pnorm for x in (-37.5, -5.657) union (5.657, 37.5) */
+                xsq = 1.0 / (x * x); /* (1./x)*(1./x) might be better */
+                xnum = p[5] * xsq;
+                xden = xsq;
+                for (i = 0; i < 4; ++i) {
+                    xnum = (xnum + p[i]) * xsq;
+                    xden = (xden + q[i]) * xsq;
+                }
+                temp = xsq * (xnum + p[4]) / (xden + q[4]);
+                temp = (jstat.ONE_SQRT_2PI - temp) / y;
+
+                /* do_del */
+                xsq = jstat.trunc(x * 16) / 16;
+                del = (x - xsq) * (x + xsq);
+                if(log_p) {
+                    cum = (-xsq * xsq * 0.5) + (-del * 0.5) + Math.log(temp);
+                    if((lower && x > 0.) || (upper && x <= 0.))
+                        ccum = jstat.log1p(-Math.exp(-xsq * xsq * 0.5) *
+                            Math.exp(-del * 0.5) * temp);
+                }
+                else {
+                    cum = Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp;
+                    ccum = 1.0 - cum;
+                }
+                /* end do_del */
+
+                /* swap_tail */
+                if (x > 0.0) {/* swap  ccum <--> cum */
+                    temp = cum;
+                    if(lower) {
+                        cum = ccum;
+
+                    }
+                    ccum = temp;
+                }
+            /* end swap_tail */
+
+            } else { /* large x such that probs are 0 or 1 */
+                if(x > 0) {
+                    cum = (log_p) ? 0.0 : 1.0;  // R_D__1
+                    ccum = (log_p) ? Number.NEGATIVE_INFINITY : 0.0;  //R_D__0;
+                } else {
+                    cum = (log_p) ? Number.NEGATIVE_INFINITY : 0.0;  //R_D__0;
+                    ccum = (log_p) ? 0.0 : 1.0;  // R_D__1
+                }
+            }
+
+            return [cum, ccum];
+        }
+
+        var p, cp;
+        var mu = this._mean;
+        var sigma = this._sigma;
+        var R_DT_0, R_DT_1;
+
+        if(lower_tail) {
+            if(log_p) {
+                R_DT_0 = Number.NEGATIVE_INFINITY;
+                R_DT_1 = 0.0;
+            } else {
+                R_DT_0 = 0.0;
+                R_DT_1 = 1.0;
+            }
+        } else {
+            if(log_p) {
+                R_DT_0 = 0.0;
+                R_DT_1 = Number.NEGATIVE_INFINITY;
+            } else {
+                R_DT_0 = 1.0;
+                R_DT_1 = 0.0;
+            }
+        }
+
+        if(!jstat.isFinite(x) && mu == x) return Number.NaN;
+        if(sigma <= 0) {
+            if(sigma < 0) {
+                console.warn("Sigma is less than 0");
+                return Number.NaN;
+            }
+            return (x < mu) ? R_DT_0 : R_DT_1;
+        }
+
+        p = (x - mu) / sigma;
+
+        if(!jstat.isFinite(p)) {
+            return (x < mu) ? R_DT_0 : R_DT_1;
+        }
+
+        x = p;
+
+        // pnorm_both(x, &p, &cp, (lower_tail ? 0 : 1), log_p);
+        // result[0] == &p
+        // result[1] == &cp
+
+        var result = pnorm_both(x, p, cp, (lower_tail ? false : true), log_p);
+
+        return (lower_tail ? result[0] : result[1]);
+
+    },
+    getMean: function() {
+        return this._mean;
+    },
+    getSigma: function() {
+        return this._sigma;
+    },
+    getVariance: function() {
+        return this._sigma*this._sigma;
+    }
+});
+
+/**
+ *  A Log-normal distribution object
+ */
+var LogNormalDistribution = ContinuousDistribution.extend({
+    init: function(location, scale) {
+        this._super('LogNormal')
+        this._location = parseFloat(location);
+        this._scale = parseFloat(scale);
+        this._string = "LogNormal ("+this._location.toFixed(2)+", " + this._scale.toFixed(2) + ")";
+    },
+    _pdf: function(x, give_log) {
+        var y;
+        var sdlog = this._scale;
+        var meanlog = this._location;
+        if(give_log == null) {
+            give_log = false;
+        }
+
+        if(sdlog <= 0) throw "Illegal parameter in _pdf";
+
+        if(x <= 0) {
+            return (give_log) ? Number.NEGATIVE_INFINITY : 0.0;
+        }
+
+        y = (Math.log(x) - meanlog) / sdlog;
+
+        return (give_log ? -(jstat.LN_SQRT_2PI + 0.5 * y * y + Math.log(x * sdlog)) :
+            jstat.ONE_SQRT_2PI * Math.exp(-0.5 * y * y) / (x * sdlog));
+
+    },
+    _cdf: function(x, lower_tail, log_p) {
+        var sdlog = this._scale;
+        var meanlog = this._location;
+        if(lower_tail == null) {
+            lower_tail = true;
+        }
+        if(log_p == null) {
+            log_p = false;
+        }
+
+
+        if(sdlog <= 0) {
+            throw "illegal std in _cdf";
+        }
+
+        if(x > 0) {
+            var nd = new NormalDistribution(meanlog, sdlog);
+            return nd._cdf(Math.log(x), lower_tail, log_p);
+        }
+        if(lower_tail) {
+            return (log_p) ? Number.NEGATIVE_INFINITY : 0.0;    // R_D__0
+        } else {
+            return (log_p) ? 0.0 : 1.0;                         // R_D__1
+        }
+    },
+    getLocation: function() {
+        return this._location;
+    },
+    getScale: function() {
+        return this._scale;
+    },
+    getMean: function() {
+        return Math.exp((this._location + this._scale) /  2);
+    },
+    getVariance: function() {
+        var ans = (Math.exp(this._scale)-1)*Math.exp(2*this._location+this._scale);
+        return ans;
+    }
+});
+
+
+/**
+ *  Gamma distribution object
+ */
+var GammaDistribution = ContinuousDistribution.extend({
+    init: function(shape, scale) {
+        this._super('Gamma');
+        this._shape = parseFloat(shape);
+        this._scale = parseFloat(scale);
+        this._string = "Gamma ("+this._shape.toFixed(2)+", " + this._scale.toFixed(2) + ")";
+    },
+    _pdf: function(x, give_log) {
+        var pr;
+        var shape = this._shape;
+        var scale = this._scale;
+        if(give_log == null) {
+            give_log = false;    // default value
+        }
+
+        if(shape < 0 || scale <= 0) {
+            throw "Illegal argument in _pdf";
+        }
+
+        if(x < 0) {
+            return (give_log) ? Number.NEGATIVE_INFINITY : 0.0; // R_D__0
+        }
+        if(shape == 0) { /* point mass at 0 */
+            return (x == 0) ? Number.POSITIVE_INFINITY : (give_log) ? Number.NEGATIVE_INFINITY : 0.0;   // R_D__0
+        }
+        if(x == 0) {
+            if(shape < 1) return Number.POSITIVE_INFINITY;
+            if(shape > 1) return (give_log) ? Number.NEGATIVE_INFINITY : 0.0; // R_D__0
+            /* else */
+            return (give_log) ? -Math.log(scale) : 1/scale;
+        }
+
+        if(shape < 1) {
+            pr = jstat.dopois_raw(shape, x/scale, give_log);
+            return give_log ? pr + Math.log(shape/x) : pr*shape/x;
+        }
+        /* else shape >= 1 */
+        pr = jstat.dopois_raw(shape-1, x/scale, give_log);
+        return give_log ? pr - Math.log(scale) : pr/scale;
+
+    },
+    /**
+     *	This function computes the distribution function for the
+     *	gamma distribution with shape parameter alph and scale parameter
+     *	scale.	This is also known as the incomplete gamma function.
+     *	See Abramowitz and Stegun (6.5.1) for example.
+     */
+    _cdf: function(x, lower_tail, log_p) {
+        /* define USE_PNORM */
+        function USE_PNORM() {
+            pn1 = Math.sqrt(alph) * 3.0 * (Math.pow(x/alph,1.0/3.0) + 1.0 / (9.0 * alph) - 1.0);
+            var norm_dist = new NormalDistribution(0.0, 1.0);
+            return norm_dist._cdf(pn1, lower_tail, log_p);
+        }
+
+        /* Defaults */
+        if(lower_tail == null) lower_tail = true;
+        if(log_p == null) log_p = false;
+        var alph = this._shape;
+        var scale = this._scale;
+        var xbig = 1.0e+8;
+        var xlarge = 1.0e+37;
+        var alphlimit = 1e5;
+        var pn1,pn2,pn3,pn4,pn5,pn6,arg,a,b,c,an,osum,sum,n,pearson;
+
+        if(alph <= 0. || scale <= 0.) {
+            console.warn('Invalid gamma params in _cdf');
+            return Number.NaN;
+        }
+
+        x/=scale;
+        if(isNaN(x)) return x;
+        if(x <= 0.0) {
+            // R_DT_0
+            if(lower_tail) {
+                // R_D__0
+                return (log_p) ? Number.NEGATIVE_INFINITY : 0.0;
+            } else {
+                // R_D__1
+                return (log_p) ? 0.0 : 1.0;
+            }
+        }
+
+        if(alph > alphlimit) {
+            return USE_PNORM();
+        }
+
+        if(x > xbig * alph) {
+            if(x > jstat.DBL_MAX * alph) {
+                // R_DT_1
+                if(lower_tail) {
+                    // R_D__1
+                    return (log_p) ? 0.0 : 1.0;
+                } else {
+                    // R_D__0
+                    return (log_p) ? Number.NEGATIVE_INFINITY : 0.0;
+                }
+            } else {
+                return USE_PNORM();
+            }
+        }
+
+        if(x <= 1.0 || x < alph) {
+            pearson = 1; /* use pearson's series expansion */
+            arg = alph * Math.log(x) - x - jstat.lgamma(alph + 1.0);
+
+            c = 1.0;
+            sum = 1.0;
+            a = alph;
+            do {
+                a += 1.0;
+                c *= x / a;
+                sum += c;
+            } while(c > jstat.DBL_EPSILON * sum);
+        } else { /* x >= max( 1, alph) */
+            pearson = 0;/* use a continued fraction expansion */
+            arg = alph * Math.log(x) - x - jstat.lgamma(alph);
+
+            a = 1. - alph;
+            b = a + x + 1.;
+            pn1 = 1.;
+            pn2 = x;
+            pn3 = x + 1.;
+            pn4 = x * b;
+            sum = pn3 / pn4;
+
+            for (n = 1; ; n++) {
+                a += 1.;/* =   n+1 -alph */
+                b += 2.;/* = 2(n+1)-alph+x */
+                an = a * n;
+                pn5 = b * pn3 - an * pn1;
+                pn6 = b * pn4 - an * pn2;
+                if (Math.abs(pn6) > 0.) {
+                    osum = sum;
+                    sum = pn5 / pn6;
+                    if (Math.abs(osum - sum) <= jstat.DBL_EPSILON * jstat.fmin2(1.0, sum))
+                        break;
+                }
+                pn1 = pn3;
+                pn2 = pn4;
+                pn3 = pn5;
+                pn4 = pn6;
+
+                if (Math.abs(pn5) >= xlarge) {
+                    pn1 /= xlarge;
+                    pn2 /= xlarge;
+                    pn3 /= xlarge;
+                    pn4 /= xlarge;
+                }
+            }
+        }
+        arg += Math.log(sum);
+        lower_tail = (lower_tail == pearson);
+
+        if (log_p && lower_tail)
+            return(arg);
+        /* else */
+        /* sum = exp(arg); and return   if(lower_tail) sum	else 1-sum : */
+
+        if(lower_tail) {
+            return Math.exp(arg);
+        } else {
+            if(log_p) {
+                // R_Log1_Exp(arg);
+                return (arg > -Math.LN2)  ? Math.log(-jstat.expm1(arg)) : jstat.log1p(-Math.exp(arg));
+            } else {
+                return -jstat.expm1(arg);
+            }
+        }
+    },
+    getShape: function() {
+        return this._shape;
+    },
+    getScale: function() {
+        return this._scale;
+    },
+    getMean: function() {
+        return this._shape * this._scale;
+    },
+    getVariance: function() {
+        return this._shape*Math.pow(this._scale,2);
+    }
+});
+
+/**
+ *  A Beta distribution object
+ */
+var BetaDistribution = ContinuousDistribution.extend({
+    init: function(alpha, beta) {
+        this._super('Beta');
+        this._alpha = parseFloat(alpha);
+        this._beta = parseFloat(beta);
+        this._string = "Beta ("+this._alpha.toFixed(2)+", " + this._beta.toFixed(2) + ")";
+    },
+    _pdf: function(x, give_log) {
+        if(give_log == null) give_log = false; // default;
+        var a = this._alpha;
+        var b = this._beta;
+        var lval;
+        if(a <= 0 || b <= 0) {
+            console.warn('Illegal arguments in _pdf');
+            return Number.NaN;
+        }
+        if(x < 0 || x > 1) {
+            // R_D__0
+            return (give_log) ? Number.NEGATIVE_INFINITY : 0.0;
+        }
+        if(x == 0) {
+            if(a > 1) {
+                // R_D__0
+                return (give_log) ? Number.NEGATIVE_INFINITY : 0.0;
+            }
+            if(a < 1) {
+                return Number.POSITIVE_INFINITY;
+            }
+            /*a == 1 */ return (give_log) ? Math.log(b) : b;    // R_D_val(b)
+        }
+        if(x == 1) {
+            if(b > 1) {
+                // R_D__0
+                return (give_log) ? Number.NEGATIVE_INFINITY : 0.0;
+            }
+            if(b < 1) {
+                return Number.POSITIVE_INFINITY;
+            }
+            /* b == 1 */ return (give_log) ? Math.log(a) : a;   // R_D_val(a)
+        }
+        if(a<=2||b<=2) {
+            lval = (a-1)*Math.log(x) + (b-1)*jstat.log1p(-x) - jstat.logBeta(a, b);
+        } else {
+            lval = Math.log(a+b-1) + jstat.dbinom_raw(a-1, a+b-2, x, 1-x, true);
+        }
+        //R_D_exp(lval)
+        return (give_log) ? lval : Math.exp(lval);
+    },
+    _cdf: function(x, lower_tail, log_p) {
+        if(lower_tail == null) lower_tail = true;
+        if(log_p == null) log_p = false;
+
+        var pin = this._alpha;
+        var qin = this._beta;
+
+        if(pin <= 0 || qin <= 0) {
+            console.warn('Invalid argument in _cdf');
+            return Number.NaN;
+        }
+
+        if(x <= 0) {
+            //R_DT_0;
+            if(lower_tail) {
+                // R_D__0
+                return (log_p) ? Number.NEGATIVE_INFINITY : 0.0;
+            } else {
+                // R_D__1
+                return (log_p) ? 0.1 : 1.0;
+            }
+        }
+
+        if(x >= 1){
+            // R_DT_1
+            if(lower_tail) {
+                // R_D__1
+                return (log_p) ? 0.1 : 1.0;
+            } else {
+                // R_D__0
+                return (log_p) ? Number.NEGATIVE_INFINITY : 0.0;
+            }
+        }
+
+        /* else */
+        return jstat.incompleteBeta(pin, qin, x);
+    },
+    getAlpha: function() {
+        return this._alpha;
+    },
+    getBeta: function() {
+        return this._beta;
+    },
+    getMean: function() {
+        return this._alpha / (this._alpha+ this._beta);
+    },
+    getVariance: function() {
+        var ans = (this._alpha * this._beta) / (Math.pow(this._alpha+this._beta,2)*(this._alpha+this._beta+1));
+        return ans;
+    }
+});
+
+var StudentTDistribution = ContinuousDistribution.extend({
+    init: function(degreesOfFreedom, mu) {
+        this._super('StudentT');
+        this._dof = parseFloat(degreesOfFreedom);
+
+        if(mu != null) {
+            this._mu = parseFloat(mu);
+            this._string = "StudentT ("+this._dof.toFixed(2)+", " + this._mu.toFixed(2)+ ")";
+        } else {
+            this._mu = 0.0;
+            this._string = "StudentT ("+this._dof.toFixed(2)+")";
+        }
+
+    },
+    _pdf: function(x, give_log) {
+        if(give_log == null) give_log = false;
+        if(this._mu == null) {
+            return this._dt(x, give_log);
+        } else {
+            var y = this._dnt(x, give_log);
+            if(y > 1){
+                console.warn('x:' + x + ', y: ' + y);
+            }
+            return y;
+        }
+    },
+    _cdf: function(x, lower_tail, give_log) {
+        if(lower_tail == null) lower_tail = true;
+        if(give_log == null) give_log = false;
+        if(this._mu == null) {
+            return this._pt(x, lower_tail, give_log);
+        } else {
+            return this._pnt(x, lower_tail, give_log);
+        }
+    },
+    _dt: function(x, give_log) {
+        var t,u;
+        var n = this._dof;
+        if (n <= 0){
+            console.warn('Invalid parameters in _dt');
+            return Number.NaN;
+        }
+        if(!jstat.isFinite(x)) {
+            return (give_log) ? Number.NEGATIVE_INFINITY : 0.0; // R_D__0;
+        }
+
+        if(!jstat.isFinite(n)) {
+            var norm = new NormalDistribution(0.0, 1.0);
+            return norm.density(x, give_log);
+        }
+
+
+        t = -jstat.bd0(n/2.0,(n+1)/2.0) + jstat.stirlerr((n+1)/2.0) - jstat.stirlerr(n/2.0);
+        if ( x*x > 0.2*n )
+            u = Math.log( 1+ x*x/n ) * n/2;
+        else
+            u = -jstat.bd0(n/2.0,(n+x*x)/2.0) + x*x/2.0;
+
+        var p1 = jstat.TWO_PI *(1+x*x/n);
+        var p2 = t-u;
+
+        return (give_log) ? -0.5*Math.log(p1) + p2 : Math.exp(p2)/Math.sqrt(p1);       // R_D_fexp(M_2PI*(1+x*x/n), t-u);
+    },
+    _dnt: function(x, give_log) {
+        if(give_log == null) give_log = false;
+        var df = this._dof;
+        var ncp = this._mu;
+        var u;
+
+        if(df <= 0.0) {
+            console.warn("Illegal arguments _dnf");
+            return Number.NaN;
+        }
+        if(ncp == 0.0) {
+            return this._dt(x, give_log);
+        }
+
+        if(!jstat.isFinite(x)) {
+            // R_D__0
+            if(give_log) {
+                return Number.NEGATIVE_INFINITY;
+            } else {
+                return 0.0;
+            }
+        }
+
+        /* If infinite df then the density is identical to a
+           normal distribution with mean = ncp.  However, the formula
+           loses a lot of accuracy around df=1e9
+         */
+        if(!isFinite(df) || df > 1e8) {
+            var dist = new NormalDistribution(ncp, 1.);
+            return dist.density(x, give_log);
+        }
+
+        /* Do calculations on log scale to stabilize */
+
+        /* Consider two cases: x ~= 0 or not */
+        if (Math.abs(x) > Math.sqrt(df * jstat.DBL_EPSILON)) {
+            var newT = new StudentTDistribution(df+2, ncp);
+            u = Math.log(df) - Math.log(Math.abs(x)) +
+            Math.log(Math.abs(newT._pnt(x*Math.sqrt((df+2)/df), true, false) -
+                this._pnt(x, true, false)));
+        /* FIXME: the above still suffers from cancellation (but not horribly) */
+        }
+        else {  /* x ~= 0 : -> same value as for  x = 0 */
+            u = jstat.lgamma((df+1)/2) - jstat.lgamma(df/2)
+            - .5*(Math.log(Math.PI) + Math.log(df) + ncp*ncp);
+        }
+
+        return (give_log ? u : Math.exp(u));
+    },
+    _pt: function(x, lower_tail, log_p) {
+        if(lower_tail == null) lower_tail = true;
+        if(log_p == null) log_p = false;
+        var val, nx;
+        var n = this._dof;
+        var DT_0, DT_1;
+
+        if(lower_tail) {
+            if(log_p) {
+                DT_0 = Number.NEGATIVE_INFINITY;
+                DT_1 = 1.;
+            } else {
+                DT_0 = 0.;
+                DT_1 = 1.;
+            }
+        } else {
+            if(log_p) {
+                // not lower_tail but log_p
+                DT_0 = 0.;
+                DT_1 = Number.NEGATIVE_INFINITY;
+            } else {
+                // not lower_tail and not log_p
+                DT_0 = 1.;
+                DT_1 = 0.;
+            }
+        }
+
+        if(n <= 0.0) {
+            console.warn("Invalid T distribution _pt");
+            return Number.NaN;
+        }
+        var norm = new NormalDistribution(0,1);
+        if(!jstat.isFinite(x)) {
+            return (x < 0) ? DT_0 : DT_1;
+        }
+        if(!jstat.isFinite(n)) {
+            return norm._cdf(x, lower_tail, log_p);
+        }
+
+        if (n > 4e5) { /*-- Fixme(?): test should depend on `n' AND `x' ! */
+            /* Approx. from	 Abramowitz & Stegun 26.7.8 (p.949) */
+            val = 1./(4.*n);
+            return norm._cdf(x*(1. - val)/sqrt(1. + x*x*2.*val), lower_tail, log_p);
+        }
+
+        nx = 1 + (x/n)*x;
+        /* FIXME: This test is probably losing rather than gaining precision,
+         * now that pbeta(*, log_p = TRUE) is much better.
+         * Note however that a version of this test *is* needed for x*x > D_MAX */
+        if(nx > 1e100) { /* <==>  x*x > 1e100 * n  */
+            /* Danger of underflow. So use Abramowitz & Stegun 26.5.4
+	   pbeta(z, a, b) ~ z^a(1-z)^b / aB(a,b) ~ z^a / aB(a,b),
+	   with z = 1/nx,  a = n/2,  b= 1/2 :
+             */
+            var lval;
+            lval = -0.5*n*(2*Math.log(Math.abs(x)) - Math.log(n))
+            - jstat.logBeta(0.5*n, 0.5) - Math.log(0.5*n);
+            val = log_p ? lval : Math.exp(lval);
+        } else {
+            /*
+            val = (n > x * x)
+            //    ? pbeta (x * x / (n + x * x), 0.5, n / 2., 0, log_p)
+           // : pbeta (1. / nx,             n / 2., 0.5, 1, log_p);
+             */
+            if(n > x * x) {
+                var beta = new BetaDistribution(0.5, n/2.);
+                return beta._cdf(x*x/ (n + x * x), false, log_p);
+            } else {
+                beta = new BetaDistribution(n / 2., 0.5);
+                return beta._cdf(1. / nx, true, log_p);
+            }
+
+
+        }
+
+        /* Use "1 - v"  if	lower_tail  and	 x > 0 (but not both):*/
+        if(x <= 0.)
+            lower_tail = !lower_tail;
+
+        if(log_p) {
+            if(lower_tail) return jstat.log1p(-0.5*Math.exp(val));
+            else return val - M_LN2; /* = log(.5* pbeta(....)) */
+        }
+        else {
+            val /= 2.;
+            if(lower_tail) {
+                return (0.5 - val + 0.5);
+            } else {
+                return val;
+            }
+        }
+    },
+    _pnt: function(t, lower_tail, log_p) {
+
+        var dof = this._dof;
+        var ncp = this._mu;
+        var DT_0, DT_1;
+
+        if(lower_tail) {
+            if(log_p) {
+                DT_0 = Number.NEGATIVE_INFINITY;
+                DT_1 = 1.;
+            } else {
+                DT_0 = 0.;
+                DT_1 = 1.;
+            }
+        } else {
+            if(log_p) {
+                // not lower_tail but log_p
+                DT_0 = 0.;
+                DT_1 = Number.NEGATIVE_INFINITY;
+            } else {
+                // not lower_tail and not log_p
+                DT_0 = 1.;
+                DT_1 = 0.;
+            }
+        }
+
+        var albeta, a, b, del, errbd, lambda, rxb, tt, x;
+        var geven, godd, p, q, s, tnc, xeven, xodd;
+        var it, negdel;
+
+        /* note - itrmax and errmax may be changed to suit one's needs. */
+        var ITRMAX = 1000;
+        var ERRMAX = 1.e-7;
+
+        if(dof <= 0.0) {
+            return Number.NaN;
+        } else if (dof == 0.0) {
+            return this._pt(t);
+        }
+
+        if(!jstat.isFinite(t)) {
+            return (t < 0) ? DT_0 : DT_1;
+        }
+        if(t >= 0.) {
+            negdel = false;
+            tt = t;
+            del = ncp;
+        } else {
+            /* We deal quickly with left tail if extreme,
+	   since pt(q, df, ncp) <= pt(0, df, ncp) = \Phi(-ncp) */
+            if(ncp >= 40 && (!log_p || !lower_tail)) {
+                return DT_0;
+            }
+            negdel = true;
+            tt = -t;
+            del = -ncp;
+        }
+
+        if(dof > 4e5 || del*del > 2* Math.LN2 * (-(jstat.DBL_MIN_EXP))) {
+            /*-- 2nd part: if del > 37.62, then p=0 below
+	    FIXME: test should depend on `df', `tt' AND `del' ! */
+            /* Approx. from	 Abramowitz & Stegun 26.7.10 (p.949) */
+            s=1./(4.*dof);
+            var norm = new NormalDistribution(del, Math.sqrt(1. + tt*tt*2.*s));
+            var result = norm._cdf(tt*(1.-s), lower_tail != negdel, log_p);
+            return result;
+        }
+
+        /* initialize twin series */
+        /* Guenther, J. (1978). Statist. Computn. Simuln. vol.6, 199. */
+        x = t * t;
+        rxb = dof/(x + dof);/* := (1 - x) {x below} -- but more accurately */
+        x = x / (x + dof);/* in [0,1) */
+        if (x > 0.) {/* <==>  t != 0 */
+            lambda = del * del;
+            p = .5 * Math.exp(-.5 * lambda);
+            if(p == 0.) {   // underflow!
+                console.warn("underflow in _pnt");
+                return DT_0;
+            }
+            q = jstat.SQRT_2dPI * p * del;
+            s = .5 - p;
+            if(s < 1e-7) {
+                s = -0.5 * jstat.expm1(-0.5 * lambda);
+            }
+            a = .5;
+            b = .5 * dof;
+            /* rxb = (1 - x) ^ b   [ ~= 1 - b*x for tiny x --> see 'xeven' below]
+             * where '(1 - x)' =: rxb {accurately!} above */
+            rxb = Math.pow(rxb, b);
+            albeta = jstat.LN_SQRT_PI + jstat.lgamma(b) - jstat.lgamma(.5 + b);
+            /* TODO: change incompleteBeta function to accept lower_tail and p_log */
+            xodd = jstat.incompleteBeta(a, b, x);
+            godd = 2. * rxb * Math.exp(a * Math.log(x) - albeta);
+            tnc = b * x;
+            xeven = (tnc < jstat.DBL_EPSILON) ? tnc : 1. - rxb;
+            geven = tnc * rxb;
+            tnc = p * xodd + q * xeven;
+
+            /* repeat until convergence or iteration limit */
+            for(it = 1; it <= ITRMAX; it++) {
+                a += 1.;
+                xodd  -= godd;
+                xeven -= geven;
+                godd  *= x * (a + b - 1.) / a;
+                geven *= x * (a + b - .5) / (a + .5);
+                p *= lambda / (2 * it);
+                q *= lambda / (2 * it + 1);
+                tnc += p * xodd + q * xeven;
+                s -= p;
+                /* R 2.4.0 added test for rounding error here. */
+                if(s < -1.e-10) { /* happens e.g. for (t,df,ncp)=(40,10,38.5), after 799 it.*/
+                    //console.write("precision error _pnt");
+                    break;
+                }
+                if(s <= 0 && it > 1) break;
+
+                errbd = 2. * s * (xodd - godd);
+
+                if(Math.abs(errbd) < ERRMAX) break;/*convergence*/
+            }
+
+            if(it == ITRMAX) {
+                throw "Non-convergence _pnt";
+            }
+        } else {
+            tnc = 0.;
+        }
+        norm = new NormalDistribution(0,1);
+        tnc += norm._cdf(-del, /*lower*/true, /*log_p*/ false);
+
+        lower_tail = lower_tail != negdel; /* xor */
+        if(tnc > 1 - 1e-10 && lower_tail) {
+            //console.warn("precision error _pnt");
+        }
+        var res = jstat.fmin2(tnc, 1.);
+        if(lower_tail) {
+            if(log_p) {
+                return Math.log(res);
+            } else {
+                return res;
+            }
+        } else {
+            if(log_p) {
+                return jstat.log1p(-(res));
+            } else {
+                return (0.5 - (res) + 0.5);
+            }
+        }
+    },
+    getDegreesOfFreedom: function() {
+        return this._dof;
+    },
+    getNonCentralityParameter: function() {
+        return this._mu;
+    },
+    getMean: function() {
+        if(this._dof > 1) {
+            var ans = (1/2)*Math.log(this._dof/2) + jstat.lgamma((this._dof-1)/2) - jstat.lgamma(this._dof/2)
+            return Math.exp(ans) * this._mu;
+        } else {
+            return Number.NaN;
+        }
+    },
+    getVariance: function() {
+        if(this._dof > 2) {
+            var ans = this._dof * (1 + this._mu*this._mu)/(this._dof-2) - (((this._mu*this._mu * this._dof) / 2) * Math.pow(Math.exp(jstat.lgamma((this._dof - 1)/2)-jstat.lgamma(this._dof/2)), 2));
+            return ans;
+        } else {
+            return Number.NaN;
+        }
+    }
+});
+
+
+/******************************************************************************/
+/*                      jQuery Flot graph objects                             */
+/******************************************************************************/
+
+var Plot = Class.extend({
+    init: function(id, options) {
+        this._container = '#' + String(id);
+        this._plots = [];
+        this._flotObj = null;
+        this._locked = false;
+
+        if(options != null) {
+            this._options = options;
+        } else {
+            this._options = {
+            };
+        }
+
+    },
+    getContainer: function() {
+        return this._container;
+    },
+    getGraph: function() {
+        return this._flotObj;
+    },
+    setData: function(data) {
+        this._plots = data;
+    },
+    clear: function() {
+        this._plots = [];
+    //this.render();
+    },
+    showLegend: function() {
+        this._options.legend = {
+            show: true
+        }
+        this.render();
+    },
+    hideLegend: function() {
+        this._options.legend = {
+            show: false
+        }
+        this.render();
+    },
+    render: function() {
+        this._flotObj = null;
+        this._flotObj = $.plot($(this._container), this._plots, this._options);
+    }
+});
+
+var DistributionPlot = Plot.extend({
+    init: function(id, distribution, range, options) {
+        this._super(id, options);
+        this._showPDF = true;
+        this._showCDF = false;
+        this._pdfValues = [];     // raw values for pdf
+        this._cdfValues = [];     // raw values for cdf
+        this._maxY = 1;
+        this._plotType = 'line';    // line, point, both
+        this._fill = false;
+        this._distribution = distribution;    // underlying PDF
+        // Range object for the plot
+        if(range != null && Range.validate(range)) {
+            this._range = range;
+        } else {
+            this._range = this._distribution.getRange(); // no range supplied, use distribution default
+        }
+
+        // render
+        if(this._distribution != null) {
+            this._maxY = this._generateValues();   // create the pdf/cdf values in the ctor
+        } else {
+            this._options.xaxis = {
+                min: range.getMinimum(),
+                max: range.getMaximum()
+            }
+            this._options.yaxis = {
+                max: 1
+            }
+        }
+
+
+
+        this.render();
+    },
+    setHover: function(bool) {
+        if(bool) {
+            if(this._options.grid == null) {
+                this._options.grid = {
+                    hoverable: true,
+                    mouseActiveRadius: 25
+                }
+            } else {
+                this._options.grid.hoverable = true,
+                this._options.grid.mouseActiveRadius = 25
+            }
+            function showTooltip(x, y, contents, color) {
+                $('<div id="jstat_tooltip">' + contents + '</div>').css( {
+                    position: 'absolute',
+                    display: 'none',
+                    top: y + 15,
+                    'font-size': 'small',
+                    left: x + 5,
+                    border: '1px solid ' + color[1],
+                    color: color[2],
+                    padding: '5px',
+                    'background-color': color[0],
+                    opacity: 0.80
+                }).appendTo("body").show();
+            }
+            var previousPoint = null;
+            $(this._container).bind("plothover", function(event, pos, item) {
+                $("#x").text(pos.x.toFixed(2));
+                $("#y").text(pos.y.toFixed(2));
+                if (item) {
+                    if (previousPoint != item.datapoint) {
+                        previousPoint = item.datapoint;
+                        $("#jstat_tooltip").remove();
+                        var x = jstat.toSigFig(item.datapoint[0],2), y = jstat.toSigFig(item.datapoint[1], 2);
+                        var text = null;
+                        var color = item.series.color;
+                        if(item.series.label == 'PDF') {
+                            text = "P(" + x + ") = " + y;
+                            color = ["#fee", "#fdd", "#C05F5F"];
+                        } else {
+                            // cdf
+                            text = "F(" + x + ") = " + y;
+                            color = ["#eef", "#ddf", "#4A4AC0"];
+                        }
+                        showTooltip(item.pageX, item.pageY, text, color);
+                    }
+                }
+                else {
+                    $("#jstat_tooltip").remove();
+                    previousPoint = null;
+                }
+            });
+            $(this._container).bind("mouseleave", function() {
+                if($('#jstat_tooltip').is(':visible')) {
+                    $('#jstat_tooltip').remove();
+                    previousPoint = null;
+                }
+            });
+        } else {
+            // unbind
+            if(this._options.grid == null) {
+                this._options.grid = {
+                    hoverable: false
+                }
+            } else {
+                this._options.grid.hoverable = false
+            }
+            $(this._container).unbind("plothover");
+        }
+
+        this.render();
+    },
+    setType: function(type) {
+        this._plotType = type;
+        var lines = {};
+        var points = {};
+        if(this._plotType == 'line') {
+            lines.show = true;
+            points.show = false;
+        } else if(this._plotType == 'points') {
+            lines.show = false;
+            points.show = true;
+        } else if(this._plotType == 'both') {
+            lines.show = true;
+            points.show = true;
+        }
+        if(this._options.series == null) {
+            this._options.series = {
+                lines: lines,
+                points: points
+            }
+        } else {
+            if(this._options.series.lines == null) {
+                this._options.series.lines = lines;
+            } else {
+                this._options.series.lines.show = lines.show;
+            }
+            if(this._options.series.points == null) {
+                this._options.series.points = points;
+            } else {
+                this._options.series.points.show = points.show;
+            }
+        }
+
+        this.render();
+    },
+    setFill: function(bool) {
+        this._fill = bool;
+        if(this._options.series == null) {
+            this._options.series = {
+                lines: {
+                    fill: bool
+                }
+            }
+        } else {
+            if(this._options.series.lines == null) {
+                this._options.series.lines = {
+                    fill: bool
+                }
+            } else {
+                this._options.series.lines.fill = bool;
+            }
+        }
+        this.render();
+    },
+    clear: function() {
+        this._super();
+        this._distribution = null;
+        this._pdfValues = [];
+        this._cdfValues = [];
+        this.render();
+    },
+    _generateValues: function() {
+        this._cdfValues = [];     // reinitialize the arrays.
+        this._pdfValues = [];
+
+        var xs = this._range.getPoints();
+
+        this._options.xaxis = {
+            min: xs[0],
+            max: xs[xs.length-1]
+        }
+        var pdfs = this._distribution.density(this._range);
+        var cdfs = this._distribution.cumulativeDensity(this._range);
+        for(var i = 0; i < xs.length; i++) {
+            if(pdfs[i] == Number.POSITIVE_INFINITY || pdfs[i] == Number.NEGATIVE_INFINITY) {
+                pdfs[i] = null;
+            }
+            if(cdfs[i] == Number.POSITIVE_INFINITY || cdfs[i] == Number.NEGATIVE_INFINITY) {
+                cdfs[i] = null;
+            }
+            this._pdfValues.push([xs[i], pdfs[i]]);
+            this._cdfValues.push([xs[i], cdfs[i]]);
+        }
+        return jstat.max(pdfs);
+    },
+    showPDF: function() {
+        this._showPDF = true;
+        this.render();
+    },
+    hidePDF: function() {
+        this._showPDF = false;
+        this.render();
+    },
+    showCDF: function() {
+        this._showCDF = true;
+        this.render();
+    },
+    hideCDF: function() {
+        this._showCDF = false;
+        this.render();
+    },
+    setDistribution: function(distribution, range) {
+        this._distribution = distribution;
+        if(range != null) {
+            this._range = range;
+        } else {
+            this._range = distribution.getRange();
+        }
+        this._maxY = this._generateValues();
+        this._options.yaxis = {
+            max: this._maxY*1.1
+        }
+
+        this.render();
+    },
+    getDistribution: function() {
+        return this._distribution;
+    },
+    getRange: function() {
+        return this._range;
+    },
+    setRange: function(range) {
+        this._range = range;
+        this._generateValues();
+        this.render();
+    },
+    render: function() {
+        if(this._distribution != null) {
+            if(this._showPDF && this._showCDF) {
+                this.setData([{
+                    yaxis: 1,
+                    data:this._pdfValues,
+                    color: 'rgb(237,194,64)',
+                    clickable: false,
+                    hoverable: true,
+                    label: "PDF"
+                }, {
+                    yaxis: 2,
+                    data:this._cdfValues,
+                    clickable: false,
+                    color: 'rgb(175,216,248)',
+                    hoverable: true,
+                    label: "CDF"
+                }]);
+                this._options.yaxis = {
+                    max: this._maxY*1.1
+                }
+            } else if(this._showPDF) {
+                this.setData([{
+                    data:this._pdfValues,
+                    hoverable: true,
+                    color: 'rgb(237,194,64)',
+                    clickable: false,
+                    label: "PDF"
+                }]);
+                this._options.yaxis = {
+                    max: this._maxY*1.1
+                }
+            } else if(this._showCDF) {
+                this.setData([{
+                    data:this._cdfValues,
+                    hoverable: true,
+                    color: 'rgb(175,216,248)',
+                    clickable: false,
+                    label: "CDF"
+                }]);
+                this._options.yaxis = {
+                    max: 1.1
+                }
+            }
+        } else {
+            this.setData([]);
+        }
+        this._super();  // Call the parent plot method
+    }
+});
+
+var DistributionFactory = {};
+DistributionFactory.build = function(json) {
+    /*
+    if(json.name == null) {
+        try{
+            json = JSON.parse(json);
+        }
+        catch(err) {
+            throw "invalid JSON";
+        }
+
+    // try to parse it
+    }*/
+
+    /*
+    if(json.name != null) {
+        var name = json.name;
+    } else {
+        throw "Malformed JSON provided to DistributionBuilder " + json;
+    }
+     */
+    if(json.NormalDistribution) {
+        if(json.NormalDistribution.mean != null && json.NormalDistribution.standardDeviation != null) {
+            return new NormalDistribution(json.NormalDistribution.mean[0], json.NormalDistribution.standardDeviation[0]);
+        } else {
+            throw "Malformed JSON provided to DistributionBuilder " + json;
+        }
+    } else if (json.LogNormalDistribution) {
+        if(json.LogNormalDistribution.location != null && json.LogNormalDistribution.scale != null) {
+            return new LogNormalDistribution(json.LogNormalDistribution.location[0], json.LogNormalDistribution.scale[0]);
+        } else {
+            throw "Malformed JSON provided to DistributionBuilder " + json;
+        }
+    } else if (json.BetaDistribution) {
+        if(json.BetaDistribution.alpha != null && json.BetaDistribution.beta != null) {
+            return new BetaDistribution(json.BetaDistribution.alpha[0], json.BetaDistribution.beta[0]);
+        } else {
+            throw "Malformed JSON provided to DistributionBuilder " + json;
+        }
+    } else if (json.GammaDistribution) {
+        if(json.GammaDistribution.shape != null && json.GammaDistribution.scale != null) {
+            return new GammaDistribution(json.GammaDistribution.shape[0], json.GammaDistribution.scale[0]);
+        } else {
+            throw "Malformed JSON provided to DistributionBuilder " + json;
+        }
+    } else if (json.StudentTDistribution) {
+        if(json.StudentTDistribution.degreesOfFreedom != null && json.StudentTDistribution.nonCentralityParameter != null) {
+            return new StudentTDistribution(json.StudentTDistribution.degreesOfFreedom[0], json.StudentTDistribution.nonCentralityParameter[0]);
+        } else if(json.StudentTDistribution.degreesOfFreedom != null) {
+            return new StudentTDistribution(json.StudentTDistribution.degreesOfFreedom[0]);
+        } else {
+            throw "Malformed JSON provided to DistributionBuilder " + json;
+        }
+    } else {
+        throw "Malformed JSON provided to DistributionBuilder " + json;
+    }
+}
+
+