1"use strict";
  2exports.__esModule = true;
  3var Errors_1 = require("../Errors");
  4/**
  5 * Returns the continued fraction for the incomplete Beta function with parameters a and b modified by Lentz's method
  6 * evaluated at x. For more information see http://jstat.github.io/special-functions.html#betacf
  7 */
  8function betacf(x, a, b) {
  9    var fpmin = 1e-30;
 10    var m = 1;
 11    var qab = a + b;
 12    var qap = a + 1;
 13    var qam = a - 1;
 14    var c = 1;
 15    var d = 1 - qab * x / qap;
 16    var m2, aa, del, h;
 17    // These q's will be used in factors that occur in the coefficients
 18    if (Math.abs(d) < fpmin)
 19        d = fpmin;
 20    d = 1 / d;
 21    h = d;
 22    for (; m <= 100; m++) {
 23        m2 = 2 * m;
 24        aa = m * (b - m) * x / ((qam + m2) * (a + m2));
 25        // One step (the even one) of the recurrence
 26        d = 1 + aa * d;
 27        if (Math.abs(d) < fpmin)
 28            d = fpmin;
 29        c = 1 + aa / c;
 30        if (Math.abs(c) < fpmin)
 31            c = fpmin;
 32        d = 1 / d;
 33        h *= d * c;
 34        aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
 35        // Next step of the recurrence (the odd one)
 36        d = 1 + aa * d;
 37        if (Math.abs(d) < fpmin)
 38            d = fpmin;
 39        c = 1 + aa / c;
 40        if (Math.abs(c) < fpmin)
 41            c = fpmin;
 42        d = 1 / d;
 43        del = d * c;
 44        h *= del;
 45        if (Math.abs(del - 1.0) < 3e-7)
 46            break;
 47    }
 48    return h;
 49}
 50exports.betacf = betacf;
 51/**
 52 * Returns the incomplete Beta function evaluated at (x,a,b). See http://jstat.github.io/special-functions.html#ibeta
 53 * for more information.
 54 */
 55function ibeta(x, a, b) {
 56    // Factors in front of the continued fraction.
 57    var bt = (x === 0 || x === 1) ? 0 :
 58        Math.exp(gammaln(a + b) - gammaln(a) -
 59            gammaln(b) + a * Math.log(x) + b *
 60            Math.log(1 - x));
 61    if (x < 0 || x > 1)
 62        // WARNING: I changed this to 0, because TS complains about doing numerical operations on boolean values.
 63        // Still safe in javascript, but not TS.
 64        return 0;
 65    if (x < (a + 1) / (a + b + 2))
 66        // Use continued fraction directly.
 67        return bt * betacf(x, a, b) / a;
 68    // else use continued fraction after making the symmetry transformation.
 69    return 1 - bt * betacf(1 - x, b, a) / b;
 70}
 71exports.ibeta = ibeta;
 72/**
 73 * Returns the inverse of the incomplete Beta function evaluated at (p,a,b). For more information see
 74 * http://jstat.github.io/special-functions.html#ibetainv
 75 */
 76function ibetainv(p, a, b) {
 77    var EPS = 1e-8;
 78    var a1 = a - 1;
 79    var b1 = b - 1;
 80    var j = 0;
 81    var lna, lnb, pp, t, u, err, x, al, h, w, afac;
 82    if (p <= 0)
 83        return 0;
 84    if (p >= 1)
 85        return 1;
 86    if (a >= 1 && b >= 1) {
 87        pp = (p < 0.5) ? p : 1 - p;
 88        t = Math.sqrt(-2 * Math.log(pp));
 89        x = (2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t;
 90        if (p < 0.5)
 91            x = -x;
 92        al = (x * x - 3) / 6;
 93        h = 2 / (1 / (2 * a - 1) + 1 / (2 * b - 1));
 94        w = (x * Math.sqrt(al + h) / h) - (1 / (2 * b - 1) - 1 / (2 * a - 1)) *
 95            (al + 5 / 6 - 2 / (3 * h));
 96        x = a / (a + b * Math.exp(2 * w));
 97    }
 98    else {
 99        lna = Math.log(a / (a + b));
100        lnb = Math.log(b / (a + b));
101        t = Math.exp(a * lna) / a;
102        u = Math.exp(b * lnb) / b;
103        w = t + u;
104        if (p < t / w)
105            x = Math.pow(a * w * p, 1 / a);
106        else
107            x = 1 - Math.pow(b * w * (1 - p), 1 / b);
108    }
109    afac = -gammaln(a) - gammaln(b) + gammaln(a + b);
110    for (; j < 10; j++) {
111        if (x === 0 || x === 1)
112            return x;
113        err = ibeta(x, a, b) - p;
114        t = Math.exp(a1 * Math.log(x) + b1 * Math.log(1 - x) + afac);
115        u = err / t;
116        x -= (t = u / (1 - 0.5 * Math.min(1, u * (a1 / x - b1 / (1 - x)))));
117        if (x <= 0)
118            x = 0.5 * (x + t);
119        if (x >= 1)
120            x = 0.5 * (x + t + 1);
121        if (Math.abs(t) < EPS * x && j > 0)
122            break;
123    }
124    return x;
125}
126exports.ibetainv = ibetainv;
127/**
128 * http://jstat.github.io/distributions.html
129 */
130function inv(x, df1, df2) {
131    return df2 / (df1 * (1 / ibetainv(x, df1 / 2, df2 / 2) - 1));
132}
133exports.inv = inv;
134/**
135 * Return the standard deviation of a vector. By defaut, the population standard deviation is returned. Passing true
136 * for the flag parameter returns the sample standard deviation. See http://jstat.github.io/vector.html#stdev for
137 * more information.
138 */
139function stdev(arr, flag) {
140    return Math.sqrt(variance(arr, flag));
141}
142exports.stdev = stdev;
143/**
144 * Return the variance of a vector. By default, the population variance is calculated. Passing true as the flag
145 * indicates computes the sample variance instead. See http://jstat.github.io/vector.html#variance for more
146 * information.
147 */
148function variance(arr, flag) {
149    if ((arr.length - (flag ? 1 : 0)) === 0) {
150        throw new Errors_1.DivZeroError("Evaluation of function CORREL caused a divide by zero error.");
151    }
152    return sumsqerr(arr) / (arr.length - (flag ? 1 : 0));
153}
154exports.variance = variance;
155/**
156 * Return the sum of a vector. See http://jstat.github.io/vector.html#sum for more information.
157 */
158function sum(arr) {
159    var sum = 0;
160    var i = arr.length;
161    while (--i >= 0) {
162        sum += arr[i];
163    }
164    return sum;
165}
166exports.sum = sum;
167/**
168 * Return the mean of a vector. See http://jstat.github.io/vector.html#mean for more information.
169 */
170function mean(arr) {
171    if (arr.length === 0) {
172        throw new Errors_1.DivZeroError("Evaluation of function CORREL caused a divide by zero error.");
173    }
174    return sum(arr) / arr.length;
175}
176exports.mean = mean;
177/**
178 * Return the sum of squared errors of prediction of a vector. See http://jstat.github.io/vector.html#sumsqerr for
179 * more information.
180 */
181function sumsqerr(arr) {
182    var m = mean(arr);
183    var sum = 0;
184    var i = arr.length;
185    var tmp;
186    while (--i >= 0) {
187        tmp = arr[i] - m;
188        sum += tmp * tmp;
189    }
190    return sum;
191}
192exports.sumsqerr = sumsqerr;
193/**
194 * Return the covariance of two vectors. See http://jstat.github.io/vector.html#covariance for more information.
195 */
196function covariance(arr1, arr2) {
197    var u = mean(arr1);
198    var v = mean(arr2);
199    var arr1Len = arr1.length;
200    var sq_dev = new Array(arr1Len);
201    for (var i = 0; i < arr1Len; i++) {
202        sq_dev[i] = (arr1[i] - u) * (arr2[i] - v);
203    }
204    if ((arr1Len - 1) === 0) {
205        throw new Errors_1.DivZeroError("Evaluation of function CORREL caused a divide by zero error.");
206    }
207    return sum(sq_dev) / (arr1Len - 1);
208}
209exports.covariance = covariance;
210/**
211 * Returns the Log-Gamma function evaluated at x. See http://jstat.github.io/special-functions.html#gammaln for more
212 * information.
213 */
214function gammaln(x) {
215    var j = 0;
216    var cof = [
217        76.18009172947146, -86.50532032941677, 24.01409824083091,
218        -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5
219    ];
220    var ser = 1.000000000190015;
221    var xx, y, tmp;
222    tmp = (y = xx = x) + 5.5;
223    tmp -= (xx + 0.5) * Math.log(tmp);
224    for (; j < 6; j++)
225        ser += cof[j] / ++y;
226    return Math.log(2.5066282746310005 * ser / xx) - tmp;
227}
228exports.gammaln = gammaln;
229/**
230 * Returns the Gamma function evaluated at x. This is sometimes called the 'complete' gamma function. See
231 * http://jstat.github.io/special-functions.html#gammafn for more information.
232 */
233function gammafn(x) {
234    var p = [-1.716185138865495, 24.76565080557592, -379.80425647094563,
235        629.3311553128184, 866.9662027904133, -31451.272968848367,
236        -36144.413418691176, 66456.14382024054
237    ];
238    var q = [-30.8402300119739, 315.35062697960416, -1015.1563674902192,
239        -3107.771671572311, 22538.118420980151, 4755.8462775278811,
240        -134659.9598649693, -115132.2596755535];
241    var fact;
242    var n = 0;
243    var xden = 0;
244    var xnum = 0;
245    var y = x;
246    var i, z, yi, res;
247    if (y <= 0) {
248        res = y % 1 + 3.6e-16;
249        if (res) {
250            fact = (!(y & 1) ? 1 : -1) * Math.PI / Math.sin(Math.PI * res);
251            y = 1 - y;
252        }
253        else {
254            return Infinity;
255        }
256    }
257    yi = y;
258    if (y < 1) {
259        z = y++;
260    }
261    else {
262        z = (y -= n = (y | 0) - 1) - 1;
263    }
264    for (i = 0; i < 8; ++i) {
265        xnum = (xnum + p[i]) * z;
266        xden = xden * z + q[i];
267    }
268    res = xnum / xden + 1;
269    if (yi < y) {
270        res /= yi;
271    }
272    else if (yi > y) {
273        for (i = 0; i < n; ++i) {
274            res *= y;
275            y++;
276        }
277    }
278    if (fact) {
279        res = fact / res;
280    }
281    return res;
282}
283exports.gammafn = gammafn;
284/**
285 * Returns the value of x in the cdf of the Gamma distribution with the parameters shape (k) and scale (theta). Notice
286 * that if using the alpha beta convention, scale = 1/beta. For more information see
287 * http://jstat.github.io/distributions.html#jStat.gamma.cdf
288 */
289function cdf(x, df1, df2) {
290    return ibeta((df1 * x) / (df1 * x + df2), df1 / 2, df2 / 2);
291}
292exports.cdf = cdf;
293/**
294 * Returns the error function evaluated at x. See also:
295 *
296 * * http://jstat.github.io/special-functions.html#erf
297 *
298 * * https://github.com/jstat/jstat/search?utf8=%E2%9C%93&q=erf&type=
299 *
300 * @param x to evaluate
301 * @returns {number} error
302 */
303function erf(x) {
304    var cof = [-1.3026537197817094, 6.4196979235649026e-1, 1.9476473204185836e-2,
305        -9.561514786808631e-3, -9.46595344482036e-4, 3.66839497852761e-4,
306        4.2523324806907e-5, -2.0278578112534e-5, -1.624290004647e-6,
307        1.303655835580e-6, 1.5626441722e-8, -8.5238095915e-8,
308        6.529054439e-9, 5.059343495e-9, -9.91364156e-10,
309        -2.27365122e-10, 9.6467911e-11, 2.394038e-12,
310        -6.886027e-12, 8.94487e-13, 3.13092e-13,
311        -1.12708e-13, 3.81e-16, 7.106e-15,
312        -1.523e-15, -9.4e-17, 1.21e-16,
313        -2.8e-17];
314    var j = cof.length - 1;
315    var isneg = false;
316    var d = 0;
317    var dd = 0;
318    var t, ty, tmp, res;
319    if (x < 0) {
320        x = -x;
321        isneg = true;
322    }
323    t = 2 / (2 + x);
324    ty = 4 * t - 2;
325    for (; j > 0; j--) {
326        tmp = d;
327        d = ty * d - dd + cof[j];
328        dd = tmp;
329    }
330    res = t * Math.exp(-x * x + 0.5 * (cof[0] + ty * d) - dd);
331    return isneg ? res - 1 : 1 - res;
332}
333exports.erf = erf;
334/**
335 * Returns the value of x in the pdf of the Gamma distribution with the parameters shape (k) and scale (theta). Notice
336 * that if using the alpha beta convention, scale = 1/beta. For more information see
337 * http://jstat.github.io/distributions.html#jStat.gamma.pdf
338 */
339function pdf(x, df1, df2) {
340    if (x < 0) {
341        return undefined;
342    }
343    return Math.sqrt((Math.pow(df1 * x, df1) * Math.pow(df2, df2)) /
344        (Math.pow(df1 * x + df2, df1 + df2))) /
345        (x * betafn(df1 / 2, df2 / 2));
346}
347exports.pdf = pdf;
348function betaln(x, y) {
349    return gammaln(x) + gammaln(y) - gammaln(x + y);
350}
351exports.betaln = betaln;
352function betafn(x, y) {
353    // ensure arguments are positive
354    if (x <= 0 || y <= 0)
355        return undefined;
356    // make sure x + y doesn't exceed the upper limit of usable values
357    return (x + y > 170) ? Math.exp(betaln(x, y)) : gammafn(x) * gammafn(y) / gammafn(x + y);
358}
359exports.betafn = betafn;
360/**
361 * Cleans a float number.
362 * @param n - number to clean
363 * @returns {number} -  clean number
364 */
365function cleanFloat(n) {
366    var power = Math.pow(10, 14);
367    return Math.round(n * power) / power;
368}
369exports.cleanFloat = cleanFloat;