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