1diff --git a/README.md b/README.md
  2index d64ec50..f8f1080 100644
  3--- a/README.md
  4+++ b/README.md
  5@@ -35,12 +35,5 @@ Contingent upon cells having formats or types for primitives.
  6 ### Scrape jsdocs for functions, put in simple index.html, doc.md files to serve up simple documentation
  7 
  8 
  9-### Pull static functions outside of formulas, declare once.
 10-
 11-
 12 ### CONVERT could offer more accurate conversions for units in the same system
 13 For example 64 tbs to a qt.
 14-
 15-
 16-### Sheet.ts and parser.js should be able to concatenate criteria and values
 17-E.g. `=COUNTIFS(A7:A24, ">6", B7:B24, "<"&DATE(1969,7,20))`
 18\ No newline at end of file
 19diff --git a/src/Formulas/Math.ts b/src/Formulas/Math.ts
 20index 3528743..190f868 100644
 21--- a/src/Formulas/Math.ts
 22+++ b/src/Formulas/Math.ts
 23@@ -920,6 +920,26 @@ var SUMX2MY2 = function (arrayX, arrayY) : number {
 24 };
 25 
 26 
 27+// Private function that will recursively generate an array of the unique primitives
 28+var _countUnique = function (values: Array<any>) : Object {
 29+  var uniques = {};
 30+  for (var i = 0; i < values.length; i++) {
 31+    if (Array.isArray(values[i])) {
 32+      // For some reasons an empty range is converted to a range with a single empty string in it.
 33+      if (values[i].length === 0) {
 34+        values[i] = [""];
 35+      }
 36+      var uniquesOfArray = _countUnique(values[i]);
 37+      for (var key in uniquesOfArray) {
 38+        uniques[key] = true;
 39+      }
 40+    } else {
 41+      uniques[Serializer.serialize(values[i])] = true;
 42+    }
 43+  }
 44+  return uniques;
 45+};
 46+
 47 /**
 48  * Counts the number of unique values in a list of specified values and ranges.
 49  * @param values The values or ranges to consider for uniqueness. Supports an arbitrary number of arguments for this
 50@@ -930,27 +950,7 @@ var SUMX2MY2 = function (arrayX, arrayY) : number {
 51 var COUNTUNIQUE = function (...values) : number {
 52   ArgsChecker.checkAtLeastLength(values, 1, "COUNTUNIQUE");
 53 
 54-  // Private function that will recursively generate an array of the unique primitives
 55-  var countUniquePrivate = function (values: Array<any>) : Object {
 56-    var uniques = {};
 57-    for (var i = 0; i < values.length; i++) {
 58-      if (Array.isArray(values[i])) {
 59-        // For some reasons an empty range is converted to a range with a single empty string in it.
 60-        if (values[i].length === 0) {
 61-          values[i] = [""];
 62-        }
 63-        var uniquesOfArray = countUniquePrivate(values[i]);
 64-        for (var key in uniquesOfArray) {
 65-          uniques[key] = true;
 66-        }
 67-      } else {
 68-        uniques[Serializer.serialize(values[i])] = true;
 69-      }
 70-    }
 71-    return uniques;
 72-  };
 73-
 74-  var uniques = countUniquePrivate(values);
 75+  var uniques = _countUnique(values);
 76   return Object.keys(uniques).length;
 77 };
 78 
 79@@ -1069,11 +1069,11 @@ export {
 80   ROUND,
 81   ROUNDDOWN,
 82   ROUNDUP,
 83-  SUMPRODUCT, // Array?
 84+  SUMPRODUCT,
 85   SUMIF,
 86   SUMSQ,
 87-  SUMX2MY2, // Array?
 88-  SUMX2PY2, // Array?
 89+  SUMX2MY2,
 90+  SUMX2PY2,
 91   FLOOR,
 92   IF,
 93   COUNTIF,
 94diff --git a/src/Formulas/Statistical.ts b/src/Formulas/Statistical.ts
 95index 7e0a894..a46657f 100644
 96--- a/src/Formulas/Statistical.ts
 97+++ b/src/Formulas/Statistical.ts
 98@@ -17,6 +17,13 @@ import {
 99   SUM,
100   ABS
101 } from "./Math";
102+import {
103+  cdf,
104+  covariance,
105+  inv,
106+  pdf,
107+  stdev,
108+} from "../Utilities/MathHelpers";
109 
110 
111 /**
112@@ -189,6 +196,7 @@ var AVERAGEA = function (...values) {
113   return result / count;
114 };
115 
116+
117 /**
118  * Calculates r, the Pearson product-moment correlation coefficient of a dataset. Any text encountered in the arguments
119  * will be ignored. CORREL is synonymous with PEARSON.
120@@ -198,80 +206,6 @@ var AVERAGEA = function (...values) {
121  * @constructor
122  */
123 var CORREL = function (dataY, dataX) : number {
124-  /**
125-   * Return the standard deviation of a vector. By defaut, the population standard deviation is returned. Passing true
126-   * for the flag parameter returns the sample standard deviation. See http://jstat.github.io/vector.html#stdev for
127-   * more information.
128-   */
129-  function stdev(arr, flag) {
130-    return Math.sqrt(variance(arr, flag));
131-  }
132-
133-  /**
134-   * Return the variance of a vector. By default, the population variance is calculated. Passing true as the flag
135-   * indicates computes the sample variance instead. See http://jstat.github.io/vector.html#variance for more
136-   * information.
137-   */
138-  function variance(arr, flag) {
139-    if ((arr.length - (flag ? 1 : 0)) === 0) {
140-      throw new DivZeroError("Evaluation of function CORREL caused a divide by zero error.");
141-    }
142-    return sumsqerr(arr) / (arr.length - (flag ? 1 : 0));
143-  }
144-
145-  /**
146-   * Return the sum of a vector. See http://jstat.github.io/vector.html#sum for more information.
147-   */
148-  function sum(arr) {
149-    var sum = 0;
150-    var i = arr.length;
151-    while (--i >= 0) {
152-      sum += arr[i];
153-    }
154-    return sum;
155-  }
156-  /**
157-   * Return the mean of a vector. See http://jstat.github.io/vector.html#mean for more information.
158-   */
159-  function mean(arr) {
160-    if (arr.length === 0) {
161-      throw new DivZeroError("Evaluation of function CORREL caused a divide by zero error.");
162-    }
163-    return sum(arr) / arr.length;
164-  }
165-
166-  /**
167-   * Return the sum of squared errors of prediction of a vector. See http://jstat.github.io/vector.html#sumsqerr for
168-   * more information.
169-   */
170-  function sumsqerr(arr) {
171-    var m = mean(arr);
172-    var sum = 0;
173-    var i = arr.length;
174-    var tmp;
175-    while (--i >= 0) {
176-      tmp = arr[i] - m;
177-      sum += tmp * tmp;
178-    }
179-    return sum;
180-  }
181-
182-  /**
183-   * Return the covariance of two vectors. See http://jstat.github.io/vector.html#covariance for more information.
184-   */
185-  function covariance(arr1, arr2) {
186-    var u = mean(arr1);
187-    var v = mean(arr2);
188-    var arr1Len = arr1.length;
189-    var sq_dev = new Array(arr1Len);
190-    for (var i = 0; i < arr1Len; i++) {
191-      sq_dev[i] = (arr1[i] - u) * (arr2[i] - v);
192-    }
193-    if ((arr1Len - 1) === 0) {
194-      throw new DivZeroError("Evaluation of function CORREL caused a divide by zero error.");
195-    }
196-    return sum(sq_dev) / (arr1Len - 1);
197-  }
198   ArgsChecker.checkLength(arguments, 2, "CORREL");
199   if (!Array.isArray(dataY)) {
200     dataY = [dataY];
201@@ -328,6 +262,8 @@ var EXPONDIST = function (x, lambda, cumulative) : number {
202   return (cumulative) ? cdf(x, lambda) : pdf(x, lambda);
203 };
204 
205+
206+
207 /**
208  * Calculates the left-tailed F probability distribution (degree of diversity) for two data sets with given input x.
209  * Alternately called Fisher-Snedecor distribution or Snecdor's F distribution.
210@@ -343,178 +279,7 @@ var EXPONDIST = function (x, lambda, cumulative) : number {
211  */
212 var FDIST$LEFTTAILED = function (x, degreesFreedom1, degreesFreedom2, cumulative) : number|undefined|boolean {
213   ArgsChecker.checkLength(arguments, 4, "FDIST$LEFTTAILED");
214-  /**
215-   * Returns the Log-Gamma function evaluated at x. See http://jstat.github.io/special-functions.html#gammaln for more
216-   * information.
217-   */
218-  function gammaln(x) {
219-    var j = 0;
220-    var cof = [
221-      76.18009172947146, -86.50532032941677, 24.01409824083091,
222-      -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5
223-    ];
224-    var ser = 1.000000000190015;
225-    var xx, y, tmp;
226-    tmp = (y = xx = x) + 5.5;
227-    tmp -= (xx + 0.5) * Math.log(tmp);
228-    for (; j < 6; j++)
229-      ser += cof[j] / ++y;
230-    return Math.log(2.5066282746310005 * ser / xx) - tmp;
231-  }
232-
233-  /**
234-   * Returns the Gamma function evaluated at x. This is sometimes called the 'complete' gamma function. See
235-   * http://jstat.github.io/special-functions.html#gammafn for more information.
236-   */
237-  function gammafn(x) {
238-    var p = [-1.716185138865495, 24.76565080557592, -379.80425647094563,
239-      629.3311553128184, 866.9662027904133, -31451.272968848367,
240-      -36144.413418691176, 66456.14382024054
241-    ];
242-    var q = [-30.8402300119739, 315.35062697960416, -1015.1563674902192,
243-      -3107.771671572311, 22538.118420980151, 4755.8462775278811,
244-      -134659.9598649693, -115132.2596755535];
245-    var fact;
246-    var n = 0;
247-    var xden = 0;
248-    var xnum = 0;
249-    var y = x;
250-    var i, z, yi, res;
251-    if (y <= 0) {
252-      res = y % 1 + 3.6e-16;
253-      if (res) {
254-        fact = (!(y & 1) ? 1 : -1) * Math.PI / Math.sin(Math.PI * res);
255-        y = 1 - y;
256-      } else {
257-        return Infinity;
258-      }
259-    }
260-    yi = y;
261-    if (y < 1) {
262-      z = y++;
263-    } else {
264-      z = (y -= n = (y | 0) - 1) - 1;
265-    }
266-    for (i = 0; i < 8; ++i) {
267-      xnum = (xnum + p[i]) * z;
268-      xden = xden * z + q[i];
269-    }
270-    res = xnum / xden + 1;
271-    if (yi < y) {
272-      res /= yi;
273-    } else if (yi > y) {
274-      for (i = 0; i < n; ++i) {
275-        res *= y;
276-        y++;
277-      }
278-    }
279-    if (fact) {
280-      res = fact / res;
281-    }
282-    return res;
283-  }
284-
285-  /**
286-   * Returns the continued fraction for the incomplete Beta function with parameters a and b modified by Lentz's method
287-   * evaluated at x. For more information see http://jstat.github.io/special-functions.html#betacf.
288-   */
289-  function betacf(x, a, b) {
290-    var fpmin = 1e-30;
291-    var m = 1;
292-    var qab = a + b;
293-    var qap = a + 1;
294-    var qam = a - 1;
295-    var c = 1;
296-    var d = 1 - qab * x / qap;
297-    var m2, aa, del, h;
298-
299-    // These q's will be used in factors that occur in the coefficients
300-    if (Math.abs(d) < fpmin)
301-      d = fpmin;
302-    d = 1 / d;
303-    h = d;
304-
305-    for (; m <= 100; m++) {
306-      m2 = 2 * m;
307-      aa = m * (b - m) * x / ((qam + m2) * (a + m2));
308-      // One step (the even one) of the recurrence
309-      d = 1 + aa * d;
310-      if (Math.abs(d) < fpmin)
311-        d = fpmin;
312-      c = 1 + aa / c;
313-      if (Math.abs(c) < fpmin)
314-        c = fpmin;
315-      d = 1 / d;
316-      h *= d * c;
317-      aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
318-      // Next step of the recurrence (the odd one)
319-      d = 1 + aa * d;
320-      if (Math.abs(d) < fpmin)
321-        d = fpmin;
322-      c = 1 + aa / c;
323-      if (Math.abs(c) < fpmin)
324-        c = fpmin;
325-      d = 1 / d;
326-      del = d * c;
327-      h *= del;
328-      if (Math.abs(del - 1.0) < 3e-7)
329-        break;
330-    }
331-
332-    return h;
333-  }
334 
335-  /**
336-   * Returns the incomplete Beta function evaluated at (x,a,b). See http://jstat.github.io/special-functions.html#ibeta
337-   * for more information.
338-   */
339-  function ibeta(x, a, b) {
340-    // Factors in front of the continued fraction.
341-    var bt = (x === 0 || x === 1) ?  0 :
342-      Math.exp(gammaln(a + b) - gammaln(a) -
343-        gammaln(b) + a * Math.log(x) + b *
344-        Math.log(1 - x));
345-    if (x < 0 || x > 1)
346-      return false;
347-    if (x < (a + 1) / (a + b + 2))
348-    // Use continued fraction directly.
349-      return bt * betacf(x, a, b) / a;
350-    // else use continued fraction after making the symmetry transformation.
351-    return 1 - bt * betacf(1 - x, b, a) / b;
352-  }
353-
354-  /**
355-   * Returns the value of x in the cdf of the Gamma distribution with the parameters shape (k) and scale (theta). Notice
356-   * that if using the alpha beta convention, scale = 1/beta. For more information see
357-   * http://jstat.github.io/distributions.html#jStat.gamma.cdf
358-   */
359-  function cdf(x, df1, df2) {
360-    return ibeta((df1 * x) / (df1 * x + df2), df1 / 2, df2 / 2);
361-  }
362-
363-  /**
364-   * Returns the value of x in the pdf of the Gamma distribution with the parameters shape (k) and scale (theta). Notice
365-   * that if using the alpha beta convention, scale = 1/beta. For more information see
366-   * http://jstat.github.io/distributions.html#jStat.gamma.pdf
367-   */
368-  function pdf(x, df1, df2) {
369-    if (x < 0) {
370-      return undefined;
371-    }
372-    return Math.sqrt((Math.pow(df1 * x, df1) * Math.pow(df2, df2)) /
373-        (Math.pow(df1 * x + df2, df1 + df2))) /
374-      (x * betafn(df1/2, df2/2));
375-  }
376-  function betaln(x, y) {
377-    return gammaln(x) + gammaln(y) - gammaln(x + y);
378-  }
379-  function betafn(x, y) {
380-    // ensure arguments are positive
381-    if (x <= 0 || y <= 0)
382-      return undefined;
383-    // make sure x + y doesn't exceed the upper limit of usable values
384-    return (x + y > 170) ? Math.exp(betaln(x, y)) : gammafn(x) * gammafn(y) / gammafn(x + y);
385-  }
386   x = TypeConverter.firstValueAsNumber(x);
387   if (x < 0) {
388     throw new NumError("Function F.DIST parameter 1 value is " + x + ". It should be greater than or equal to 0.");
389@@ -525,6 +290,7 @@ var FDIST$LEFTTAILED = function (x, degreesFreedom1, degreesFreedom2, cumulative
390   return (cum) ? cdf(x, d1, d2) : pdf(x, d1, d2);
391 };
392 
393+
394 /**
395  * Returns the inverse of the (right-tailed) F probability distribution. If p = FDIST(x,...), then FINV(p,...) = x. The
396  * F distribution can be used in an F-test that compares the degree of variability in two data sets.
397@@ -536,156 +302,7 @@ var FDIST$LEFTTAILED = function (x, degreesFreedom1, degreesFreedom2, cumulative
398  */
399 var FINV = function (probability, degFreedom1, degFreedom2) : number {
400   ArgsChecker.checkLength(arguments, 3, "FINV");
401-  /**
402-   * Returns the continued fraction for the incomplete Beta function with parameters a and b modified by Lentz's method
403-   * evaluated at x. For more information see http://jstat.github.io/special-functions.html#betacf
404-   */
405-  function betacf(x, a, b) {
406-    var fpmin = 1e-30;
407-    var m = 1;
408-    var qab = a + b;
409-    var qap = a + 1;
410-    var qam = a - 1;
411-    var c = 1;
412-    var d = 1 - qab * x / qap;
413-    var m2, aa, del, h;
414-
415-    // These q's will be used in factors that occur in the coefficients
416-    if (Math.abs(d) < fpmin)
417-      d = fpmin;
418-    d = 1 / d;
419-    h = d;
420-
421-    for (; m <= 100; m++) {
422-      m2 = 2 * m;
423-      aa = m * (b - m) * x / ((qam + m2) * (a + m2));
424-      // One step (the even one) of the recurrence
425-      d = 1 + aa * d;
426-      if (Math.abs(d) < fpmin)
427-        d = fpmin;
428-      c = 1 + aa / c;
429-      if (Math.abs(c) < fpmin)
430-        c = fpmin;
431-      d = 1 / d;
432-      h *= d * c;
433-      aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
434-      // Next step of the recurrence (the odd one)
435-      d = 1 + aa * d;
436-      if (Math.abs(d) < fpmin)
437-        d = fpmin;
438-      c = 1 + aa / c;
439-      if (Math.abs(c) < fpmin)
440-        c = fpmin;
441-      d = 1 / d;
442-      del = d * c;
443-      h *= del;
444-      if (Math.abs(del - 1.0) < 3e-7)
445-        break;
446-    }
447-
448-    return h;
449-  }
450-
451-  /**
452-   * Returns the incomplete Beta function evaluated at (x,a,b). See http://jstat.github.io/special-functions.html#ibeta
453-   * for more information.
454-   */
455-  function ibeta(x, a, b) : number {
456-    // Factors in front of the continued fraction.
457-    var bt = (x === 0 || x === 1) ?  0 :
458-      Math.exp(gammaln(a + b) - gammaln(a) -
459-        gammaln(b) + a * Math.log(x) + b *
460-        Math.log(1 - x));
461-    if (x < 0 || x > 1)
462-    // WARNING: I changed this to 0, because TS complains about doing numerical operations on boolean values.
463-    // Still safe in javascript, but not TS.
464-      return 0;
465-    if (x < (a + 1) / (a + b + 2))
466-    // Use continued fraction directly.
467-      return bt * betacf(x, a, b) / a;
468-    // else use continued fraction after making the symmetry transformation.
469-    return 1 - bt * betacf(1 - x, b, a) / b;
470-  }
471-
472-  /**
473-   * Returns the Log-Gamma function evaluated at x. For more information see
474-   * http://jstat.github.io/special-functions.html#gammaln
475-   */
476-  function gammaln(x) {
477-    var j = 0;
478-    var cof = [
479-      76.18009172947146, -86.50532032941677, 24.01409824083091,
480-      -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5
481-    ];
482-    var ser = 1.000000000190015;
483-    var xx, y, tmp;
484-    tmp = (y = xx = x) + 5.5;
485-    tmp -= (xx + 0.5) * Math.log(tmp);
486-    for (; j < 6; j++)
487-      ser += cof[j] / ++y;
488-    return Math.log(2.5066282746310005 * ser / xx) - tmp;
489-  }
490 
491-  /**
492-   * Returns the inverse of the incomplete Beta function evaluated at (p,a,b). For more information see
493-   * http://jstat.github.io/special-functions.html#ibetainv
494-   */
495-  function ibetainv(p, a, b) {
496-    var EPS = 1e-8;
497-    var a1 = a - 1;
498-    var b1 = b - 1;
499-    var j = 0;
500-    var lna, lnb, pp, t, u, err, x, al, h, w, afac;
501-    if (p <= 0)
502-      return 0;
503-    if (p >= 1)
504-      return 1;
505-    if (a >= 1 && b >= 1) {
506-      pp = (p < 0.5) ? p : 1 - p;
507-      t = Math.sqrt(-2 * Math.log(pp));
508-      x = (2.30753 + t * 0.27061) / (1 + t* (0.99229 + t * 0.04481)) - t;
509-      if (p < 0.5)
510-        x = -x;
511-      al = (x * x - 3) / 6;
512-      h = 2 / (1 / (2 * a - 1)  + 1 / (2 * b - 1));
513-      w = (x * Math.sqrt(al + h) / h) - (1 / (2 * b - 1) - 1 / (2 * a - 1)) *
514-        (al + 5 / 6 - 2 / (3 * h));
515-      x = a / (a + b * Math.exp(2 * w));
516-    } else {
517-      lna = Math.log(a / (a + b));
518-      lnb = Math.log(b / (a + b));
519-      t = Math.exp(a * lna) / a;
520-      u = Math.exp(b * lnb) / b;
521-      w = t + u;
522-      if (p < t / w)
523-        x = Math.pow(a * w * p, 1 / a);
524-      else
525-        x = 1 - Math.pow(b * w * (1 - p), 1 / b);
526-    }
527-    afac = -gammaln(a) - gammaln(b) + gammaln(a + b);
528-    for(; j < 10; j++) {
529-      if (x === 0 || x === 1)
530-        return x;
531-      err = ibeta(x, a, b) - p;
532-      t = Math.exp(a1 * Math.log(x) + b1 * Math.log(1 - x) + afac);
533-      u = err / t;
534-      x -= (t = u / (1 - 0.5 * Math.min(1, u * (a1 / x - b1 / (1 - x)))));
535-      if (x <= 0)
536-        x = 0.5 * (x + t);
537-      if (x >= 1)
538-        x = 0.5 * (x + t + 1);
539-      if (Math.abs(t) < EPS * x && j > 0)
540-        break;
541-    }
542-    return x;
543-  }
544-
545-  /**
546-   * http://jstat.github.io/distributions.html
547-   */
548-  function inv(x, df1, df2) {
549-    return df2 / (df1 * (1 / ibetainv(x, df1 / 2, df2 / 2) - 1));
550-  }
551   probability = TypeConverter.firstValueAsNumber(probability);
552   if (probability <= 0.0 || probability > 1.0) {
553     throw new NumError("Function FINV parameter 1 value is " + probability
554diff --git a/src/Utilities/MathHelpers.ts b/src/Utilities/MathHelpers.ts
555new file mode 100644
556index 0000000..53c440a
557--- /dev/null
558+++ b/src/Utilities/MathHelpers.ts
559@@ -0,0 +1,336 @@
560+import {
561+  DivZeroError
562+} from "../Errors";
563+/**
564+ * Returns the continued fraction for the incomplete Beta function with parameters a and b modified by Lentz's method
565+ * evaluated at x. For more information see http://jstat.github.io/special-functions.html#betacf
566+ */
567+function betacf(x, a, b) {
568+  var fpmin = 1e-30;
569+  var m = 1;
570+  var qab = a + b;
571+  var qap = a + 1;
572+  var qam = a - 1;
573+  var c = 1;
574+  var d = 1 - qab * x / qap;
575+  var m2, aa, del, h;
576+
577+  // These q's will be used in factors that occur in the coefficients
578+  if (Math.abs(d) < fpmin)
579+    d = fpmin;
580+  d = 1 / d;
581+  h = d;
582+
583+  for (; m <= 100; m++) {
584+    m2 = 2 * m;
585+    aa = m * (b - m) * x / ((qam + m2) * (a + m2));
586+    // One step (the even one) of the recurrence
587+    d = 1 + aa * d;
588+    if (Math.abs(d) < fpmin)
589+      d = fpmin;
590+    c = 1 + aa / c;
591+    if (Math.abs(c) < fpmin)
592+      c = fpmin;
593+    d = 1 / d;
594+    h *= d * c;
595+    aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
596+    // Next step of the recurrence (the odd one)
597+    d = 1 + aa * d;
598+    if (Math.abs(d) < fpmin)
599+      d = fpmin;
600+    c = 1 + aa / c;
601+    if (Math.abs(c) < fpmin)
602+      c = fpmin;
603+    d = 1 / d;
604+    del = d * c;
605+    h *= del;
606+    if (Math.abs(del - 1.0) < 3e-7)
607+      break;
608+  }
609+
610+  return h;
611+}
612+
613+/**
614+ * Returns the incomplete Beta function evaluated at (x,a,b). See http://jstat.github.io/special-functions.html#ibeta
615+ * for more information.
616+ */
617+function ibeta(x, a, b) : number {
618+  // Factors in front of the continued fraction.
619+  var bt = (x === 0 || x === 1) ?  0 :
620+    Math.exp(gammaln(a + b) - gammaln(a) -
621+      gammaln(b) + a * Math.log(x) + b *
622+      Math.log(1 - x));
623+  if (x < 0 || x > 1)
624+  // WARNING: I changed this to 0, because TS complains about doing numerical operations on boolean values.
625+  // Still safe in javascript, but not TS.
626+    return 0;
627+  if (x < (a + 1) / (a + b + 2))
628+  // Use continued fraction directly.
629+    return bt * betacf(x, a, b) / a;
630+  // else use continued fraction after making the symmetry transformation.
631+  return 1 - bt * betacf(1 - x, b, a) / b;
632+}
633+
634+/**
635+ * Returns the inverse of the incomplete Beta function evaluated at (p,a,b). For more information see
636+ * http://jstat.github.io/special-functions.html#ibetainv
637+ */
638+function ibetainv(p, a, b) {
639+  var EPS = 1e-8;
640+  var a1 = a - 1;
641+  var b1 = b - 1;
642+  var j = 0;
643+  var lna, lnb, pp, t, u, err, x, al, h, w, afac;
644+  if (p <= 0)
645+    return 0;
646+  if (p >= 1)
647+    return 1;
648+  if (a >= 1 && b >= 1) {
649+    pp = (p < 0.5) ? p : 1 - p;
650+    t = Math.sqrt(-2 * Math.log(pp));
651+    x = (2.30753 + t * 0.27061) / (1 + t* (0.99229 + t * 0.04481)) - t;
652+    if (p < 0.5)
653+      x = -x;
654+    al = (x * x - 3) / 6;
655+    h = 2 / (1 / (2 * a - 1)  + 1 / (2 * b - 1));
656+    w = (x * Math.sqrt(al + h) / h) - (1 / (2 * b - 1) - 1 / (2 * a - 1)) *
657+      (al + 5 / 6 - 2 / (3 * h));
658+    x = a / (a + b * Math.exp(2 * w));
659+  } else {
660+    lna = Math.log(a / (a + b));
661+    lnb = Math.log(b / (a + b));
662+    t = Math.exp(a * lna) / a;
663+    u = Math.exp(b * lnb) / b;
664+    w = t + u;
665+    if (p < t / w)
666+      x = Math.pow(a * w * p, 1 / a);
667+    else
668+      x = 1 - Math.pow(b * w * (1 - p), 1 / b);
669+  }
670+  afac = -gammaln(a) - gammaln(b) + gammaln(a + b);
671+  for(; j < 10; j++) {
672+    if (x === 0 || x === 1)
673+      return x;
674+    err = ibeta(x, a, b) - p;
675+    t = Math.exp(a1 * Math.log(x) + b1 * Math.log(1 - x) + afac);
676+    u = err / t;
677+    x -= (t = u / (1 - 0.5 * Math.min(1, u * (a1 / x - b1 / (1 - x)))));
678+    if (x <= 0)
679+      x = 0.5 * (x + t);
680+    if (x >= 1)
681+      x = 0.5 * (x + t + 1);
682+    if (Math.abs(t) < EPS * x && j > 0)
683+      break;
684+  }
685+  return x;
686+}
687+
688+/**
689+ * http://jstat.github.io/distributions.html
690+ */
691+function inv(x, df1, df2) {
692+  return df2 / (df1 * (1 / ibetainv(x, df1 / 2, df2 / 2) - 1));
693+}
694+
695+
696+/**
697+ * Return the standard deviation of a vector. By defaut, the population standard deviation is returned. Passing true
698+ * for the flag parameter returns the sample standard deviation. See http://jstat.github.io/vector.html#stdev for
699+ * more information.
700+ */
701+function stdev(arr, flag) {
702+  return Math.sqrt(variance(arr, flag));
703+}
704+
705+/**
706+ * Return the variance of a vector. By default, the population variance is calculated. Passing true as the flag
707+ * indicates computes the sample variance instead. See http://jstat.github.io/vector.html#variance for more
708+ * information.
709+ */
710+function variance(arr, flag) {
711+  if ((arr.length - (flag ? 1 : 0)) === 0) {
712+    throw new DivZeroError("Evaluation of function CORREL caused a divide by zero error.");
713+  }
714+  return sumsqerr(arr) / (arr.length - (flag ? 1 : 0));
715+}
716+
717+/**
718+ * Return the sum of a vector. See http://jstat.github.io/vector.html#sum for more information.
719+ */
720+function sum(arr) {
721+  var sum = 0;
722+  var i = arr.length;
723+  while (--i >= 0) {
724+    sum += arr[i];
725+  }
726+  return sum;
727+}
728+/**
729+ * Return the mean of a vector. See http://jstat.github.io/vector.html#mean for more information.
730+ */
731+function mean(arr) {
732+  if (arr.length === 0) {
733+    throw new DivZeroError("Evaluation of function CORREL caused a divide by zero error.");
734+  }
735+  return sum(arr) / arr.length;
736+}
737+
738+/**
739+ * Return the sum of squared errors of prediction of a vector. See http://jstat.github.io/vector.html#sumsqerr for
740+ * more information.
741+ */
742+function sumsqerr(arr) {
743+  var m = mean(arr);
744+  var sum = 0;
745+  var i = arr.length;
746+  var tmp;
747+  while (--i >= 0) {
748+    tmp = arr[i] - m;
749+    sum += tmp * tmp;
750+  }
751+  return sum;
752+}
753+
754+/**
755+ * Return the covariance of two vectors. See http://jstat.github.io/vector.html#covariance for more information.
756+ */
757+function covariance(arr1, arr2) {
758+  var u = mean(arr1);
759+  var v = mean(arr2);
760+  var arr1Len = arr1.length;
761+  var sq_dev = new Array(arr1Len);
762+  for (var i = 0; i < arr1Len; i++) {
763+    sq_dev[i] = (arr1[i] - u) * (arr2[i] - v);
764+  }
765+  if ((arr1Len - 1) === 0) {
766+    throw new DivZeroError("Evaluation of function CORREL caused a divide by zero error.");
767+  }
768+  return sum(sq_dev) / (arr1Len - 1);
769+}
770+
771+/**
772+ * Returns the Log-Gamma function evaluated at x. See http://jstat.github.io/special-functions.html#gammaln for more
773+ * information.
774+ */
775+function gammaln(x) {
776+  var j = 0;
777+  var cof = [
778+    76.18009172947146, -86.50532032941677, 24.01409824083091,
779+    -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5
780+  ];
781+  var ser = 1.000000000190015;
782+  var xx, y, tmp;
783+  tmp = (y = xx = x) + 5.5;
784+  tmp -= (xx + 0.5) * Math.log(tmp);
785+  for (; j < 6; j++)
786+    ser += cof[j] / ++y;
787+  return Math.log(2.5066282746310005 * ser / xx) - tmp;
788+}
789+
790+/**
791+ * Returns the Gamma function evaluated at x. This is sometimes called the 'complete' gamma function. See
792+ * http://jstat.github.io/special-functions.html#gammafn for more information.
793+ */
794+function gammafn(x) {
795+  var p = [-1.716185138865495, 24.76565080557592, -379.80425647094563,
796+    629.3311553128184, 866.9662027904133, -31451.272968848367,
797+    -36144.413418691176, 66456.14382024054
798+  ];
799+  var q = [-30.8402300119739, 315.35062697960416, -1015.1563674902192,
800+    -3107.771671572311, 22538.118420980151, 4755.8462775278811,
801+    -134659.9598649693, -115132.2596755535];
802+  var fact;
803+  var n = 0;
804+  var xden = 0;
805+  var xnum = 0;
806+  var y = x;
807+  var i, z, yi, res;
808+  if (y <= 0) {
809+    res = y % 1 + 3.6e-16;
810+    if (res) {
811+      fact = (!(y & 1) ? 1 : -1) * Math.PI / Math.sin(Math.PI * res);
812+      y = 1 - y;
813+    } else {
814+      return Infinity;
815+    }
816+  }
817+  yi = y;
818+  if (y < 1) {
819+    z = y++;
820+  } else {
821+    z = (y -= n = (y | 0) - 1) - 1;
822+  }
823+  for (i = 0; i < 8; ++i) {
824+    xnum = (xnum + p[i]) * z;
825+    xden = xden * z + q[i];
826+  }
827+  res = xnum / xden + 1;
828+  if (yi < y) {
829+    res /= yi;
830+  } else if (yi > y) {
831+    for (i = 0; i < n; ++i) {
832+      res *= y;
833+      y++;
834+    }
835+  }
836+  if (fact) {
837+    res = fact / res;
838+  }
839+  return res;
840+}
841+
842+
843+/**
844+ * Returns the value of x in the cdf of the Gamma distribution with the parameters shape (k) and scale (theta). Notice
845+ * that if using the alpha beta convention, scale = 1/beta. For more information see
846+ * http://jstat.github.io/distributions.html#jStat.gamma.cdf
847+ */
848+function cdf(x, df1, df2) {
849+  return ibeta((df1 * x) / (df1 * x + df2), df1 / 2, df2 / 2);
850+}
851+
852+/**
853+ * Returns the value of x in the pdf of the Gamma distribution with the parameters shape (k) and scale (theta). Notice
854+ * that if using the alpha beta convention, scale = 1/beta. For more information see
855+ * http://jstat.github.io/distributions.html#jStat.gamma.pdf
856+ */
857+function pdf(x, df1, df2) {
858+  if (x < 0) {
859+    return undefined;
860+  }
861+  return Math.sqrt((Math.pow(df1 * x, df1) * Math.pow(df2, df2)) /
862+      (Math.pow(df1 * x + df2, df1 + df2))) /
863+    (x * betafn(df1/2, df2/2));
864+}
865+
866+function betaln(x, y) {
867+  return gammaln(x) + gammaln(y) - gammaln(x + y);
868+}
869+
870+function betafn(x, y) {
871+  // ensure arguments are positive
872+  if (x <= 0 || y <= 0)
873+    return undefined;
874+  // make sure x + y doesn't exceed the upper limit of usable values
875+  return (x + y > 170) ? Math.exp(betaln(x, y)) : gammafn(x) * gammafn(y) / gammafn(x + y);
876+}
877+
878+export {
879+  betacf,
880+  betafn,
881+  betaln,
882+  cdf,
883+  covariance,
884+  gammafn,
885+  gammaln,
886+  ibeta,
887+  ibetainv,
888+  inv,
889+  mean,
890+  pdf,
891+  stdev,
892+  sum,
893+  sumsqerr,
894+  variance
895+}
896\ No newline at end of file