1 /*
  2     Copyright 2008-2024
  3         Matthias Ehmann,
  4         Michael Gerhaeuser,
  5         Carsten Miller,
  6         Bianca Valentin,
  7         Alfred Wassermann,
  8         Peter Wilfahrt
  9 
 10     This file is part of JSXGraph.
 11 
 12     JSXGraph is free software dual licensed under the GNU LGPL or MIT License.
 13 
 14     You can redistribute it and/or modify it under the terms of the
 15 
 16       * GNU Lesser General Public License as published by
 17         the Free Software Foundation, either version 3 of the License, or
 18         (at your option) any later version
 19       OR
 20       * MIT License: https://github.com/jsxgraph/jsxgraph/blob/master/LICENSE.MIT
 21 
 22     JSXGraph is distributed in the hope that it will be useful,
 23     but WITHOUT ANY WARRANTY; without even the implied warranty of
 24     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 25     GNU Lesser General Public License for more details.
 26 
 27     You should have received a copy of the GNU Lesser General Public License and
 28     the MIT License along with JSXGraph. If not, see <https://www.gnu.org/licenses/>
 29     and <https://opensource.org/licenses/MIT/>.
 30  */
 31 
 32 /*global JXG: true, define: true*/
 33 /*jslint nomen: true, plusplus: true*/
 34 
 35 import JXG from "../jxg.js";
 36 import Mat from "./math.js";
 37 import Type from "../utils/type.js";
 38 
 39 /**
 40  * Functions for mathematical statistics. Most functions are like in the statistics package R.
 41  * @name JXG.Math.Statistics
 42  * @exports Mat.Statistics as JXG.Math.Statistics
 43  * @namespace
 44  */
 45 Mat.Statistics = {
 46     /**
 47      * Sums up all elements of the given array.
 48      * @param {Array} arr An array of numbers.
 49      * @returns {Number}
 50      * @memberof JXG.Math.Statistics
 51      */
 52     sum: function (arr) {
 53         var i,
 54             len = arr.length,
 55             res = 0;
 56 
 57         for (i = 0; i < len; i++) {
 58             res += arr[i];
 59         }
 60         return res;
 61     },
 62 
 63     /**
 64      * Multiplies all elements of the given array.
 65      * @param {Array} arr An array of numbers.
 66      * @returns {Number}
 67      * @memberof JXG.Math.Statistics
 68      */
 69     prod: function (arr) {
 70         var i,
 71             len = arr.length,
 72             res = 1;
 73 
 74         for (i = 0; i < len; i++) {
 75             res *= arr[i];
 76         }
 77         return res;
 78     },
 79 
 80     /**
 81      * Determines the mean value of the values given in an array.
 82      * @param {Array} arr
 83      * @returns {Number}
 84      * @memberof JXG.Math.Statistics
 85      */
 86     mean: function (arr) {
 87         if (arr.length > 0) {
 88             return this.sum(arr) / arr.length;
 89         }
 90 
 91         return 0.0;
 92     },
 93 
 94     /**
 95      * The median of a finite set of values is the value that divides the set
 96      * into two equal sized subsets.
 97      * @param {Array} arr The set of values.
 98      * @returns {Number}
 99      * @memberof JXG.Math.Statistics
100      */
101     median: function (arr) {
102         var tmp, len;
103 
104         if (arr.length > 0) {
105             if (ArrayBuffer.isView(arr)) {
106                 tmp = new Float64Array(arr);
107                 tmp.sort();
108             } else {
109                 tmp = arr.slice(0);
110                 tmp.sort(function (a, b) {
111                     return a - b;
112                 });
113             }
114             len = tmp.length;
115 
116             if (len & 1) {
117                 // odd
118                 return tmp[parseInt(len * 0.5, 10)];
119             }
120 
121             return (tmp[len * 0.5 - 1] + tmp[len * 0.5]) * 0.5;
122         }
123 
124         return 0.0;
125     },
126 
127     /**
128      * The P-th percentile ( 0 < P ≤ 100 ) of a list of N ordered values (sorted from least to greatest)
129      * is the smallest value in the list such that no more than P percent of the data is strictly less
130      * than the value and at least P percent of the data is less than or equal to that value. See {@link https://en.wikipedia.org/wiki/Percentile}.
131      *
132      * Here, the <i>linear interpolation between closest ranks</i> method is used.
133      * @param {Array} arr The set of values, need not be ordered.
134      * @param {Number|Array} percentile One or several percentiles
135      * @returns {Number|Array} Depending if a number or an array is the input for percentile, a number or an array containing the percentils
136      * is returned.
137      */
138     percentile: function (arr, percentile) {
139         var tmp,
140             len,
141             i,
142             p,
143             res = [],
144             per;
145 
146         if (arr.length > 0) {
147             if (ArrayBuffer.isView(arr)) {
148                 tmp = new Float64Array(arr);
149                 tmp.sort();
150             } else {
151                 tmp = arr.slice(0);
152                 tmp.sort(function (a, b) {
153                     return a - b;
154                 });
155             }
156             len = tmp.length;
157 
158             if (Type.isArray(percentile)) {
159                 p = percentile;
160             } else {
161                 p = [percentile];
162             }
163 
164             for (i = 0; i < p.length; i++) {
165                 per = len * p[i] * 0.01;
166                 if (parseInt(per, 10) === per) {
167                     res.push((tmp[per - 1] + tmp[per]) * 0.5);
168                 } else {
169                     res.push(tmp[parseInt(per, 10)]);
170                 }
171             }
172 
173             if (Type.isArray(percentile)) {
174                 return res;
175             } else {
176                 return res[0];
177             }
178         }
179 
180         return 0.0;
181     },
182 
183     /**
184      * Bias-corrected sample variance. A variance is a measure of how far a
185      * set of numbers are spread out from each other.
186      * @param {Array} arr
187      * @returns {Number}
188      * @memberof JXG.Math.Statistics
189      */
190     variance: function (arr) {
191         var m,
192             res,
193             i,
194             len = arr.length;
195 
196         if (len > 1) {
197             m = this.mean(arr);
198             res = 0;
199             for (i = 0; i < len; i++) {
200                 res += (arr[i] - m) * (arr[i] - m);
201             }
202             return res / (arr.length - 1);
203         }
204 
205         return 0.0;
206     },
207 
208     /**
209      * Determines the <strong>s</strong>tandard <strong>d</strong>eviation which shows how much
210      * variation there is from the average value of a set of numbers.
211      * @param {Array} arr
212      * @returns {Number}
213      * @memberof JXG.Math.Statistics
214      */
215     sd: function (arr) {
216         return Math.sqrt(this.variance(arr));
217     },
218 
219     /**
220      * Weighted mean value is basically the same as {@link JXG.Math.Statistics.mean} but here the values
221      * are weighted, i.e. multiplied with another value called <em>weight</em>. The weight values are given
222      * as a second array with the same length as the value array..
223      * @throws {Error} If the dimensions of the arrays don't match.
224      * @param {Array} arr Set of alues.
225      * @param {Array} w Weight values.
226      * @returns {Number}
227      * @memberof JXG.Math.Statistics
228      */
229     weightedMean: function (arr, w) {
230         if (arr.length !== w.length) {
231             throw new Error(
232                 "JSXGraph error (Math.Statistics.weightedMean): Array dimension mismatch."
233             );
234         }
235 
236         if (arr.length > 0) {
237             return this.mean(this.multiply(arr, w));
238         }
239 
240         return 0.0;
241     },
242 
243     /**
244      * Extracts the maximum value from the array.
245      * @param {Array} arr
246      * @returns {Number} The highest number from the array. It returns <tt>NaN</tt> if not every element could be
247      * interpreted as a number and <tt>-Infinity</tt> if an empty array is given or no element could be interpreted
248      * as a number.
249      * @memberof JXG.Math.Statistics
250      */
251     max: function (arr) {
252         return Math.max.apply(this, arr);
253     },
254 
255     /**
256      * Extracts the minimum value from the array.
257      * @param {Array} arr
258      * @returns {Number} The lowest number from the array. It returns <tt>NaN</tt> if not every element could be
259      * interpreted as a number and <tt>Infinity</tt> if an empty array is given or no element could be interpreted
260      * as a number.
261      * @memberof JXG.Math.Statistics
262      */
263     min: function (arr) {
264         return Math.min.apply(this, arr);
265     },
266 
267     /**
268      * Determines the lowest and the highest value from the given array.
269      * @param {Array} arr
270      * @returns {Array} The minimum value as the first and the maximum value as the second value.
271      * @memberof JXG.Math.Statistics
272      */
273     range: function (arr) {
274         return [this.min(arr), this.max(arr)];
275     },
276 
277     /**
278      * Determines the absolute value of every given value.
279      * @param {Array|Number} arr
280      * @returns {Array|Number}
281      * @memberof JXG.Math.Statistics
282      */
283     abs: function (arr) {
284         var i, len, res;
285 
286         if (Type.isArray(arr)) {
287             if (arr.map) {
288                 res = arr.map(Math.abs);
289             } else {
290                 len = arr.length;
291                 res = [];
292 
293                 for (i = 0; i < len; i++) {
294                     res[i] = Math.abs(arr[i]);
295                 }
296             }
297         } else if (ArrayBuffer.isView(arr)) {
298             res = arr.map(Math.abs);
299         } else {
300             res = Math.abs(arr);
301         }
302         return res;
303     },
304 
305     /**
306      * Adds up two (sequences of) values. If one value is an array and the other one is a number the number
307      * is added to every element of the array. If two arrays are given and the lengths don't match the shortest
308      * length is taken.
309      * @param {Array|Number} arr1
310      * @param {Array|Number} arr2
311      * @returns {Array|Number}
312      * @memberof JXG.Math.Statistics
313      */
314     add: function (arr1, arr2) {
315         var i,
316             len,
317             res = [];
318 
319         arr1 = Type.evalSlider(arr1);
320         arr2 = Type.evalSlider(arr2);
321 
322         if (Type.isArray(arr1) && Type.isNumber(arr2)) {
323             len = arr1.length;
324 
325             for (i = 0; i < len; i++) {
326                 res[i] = arr1[i] + arr2;
327             }
328         } else if (Type.isNumber(arr1) && Type.isArray(arr2)) {
329             len = arr2.length;
330 
331             for (i = 0; i < len; i++) {
332                 res[i] = arr1 + arr2[i];
333             }
334         } else if (Type.isArray(arr1) && Type.isArray(arr2)) {
335             len = Math.min(arr1.length, arr2.length);
336 
337             for (i = 0; i < len; i++) {
338                 res[i] = arr1[i] + arr2[i];
339             }
340         } else {
341             res = arr1 + arr2;
342         }
343 
344         return res;
345     },
346 
347     /**
348      * Divides two (sequences of) values. If two arrays are given and the lengths don't match the shortest length
349      * is taken.
350      * @param {Array|Number} arr1 Dividend
351      * @param {Array|Number} arr2 Divisor
352      * @returns {Array|Number}
353      * @memberof JXG.Math.Statistics
354      */
355     div: function (arr1, arr2) {
356         var i,
357             len,
358             res = [];
359 
360         arr1 = Type.evalSlider(arr1);
361         arr2 = Type.evalSlider(arr2);
362 
363         if (Type.isArray(arr1) && Type.isNumber(arr2)) {
364             len = arr1.length;
365 
366             for (i = 0; i < len; i++) {
367                 res[i] = arr1[i] / arr2;
368             }
369         } else if (Type.isNumber(arr1) && Type.isArray(arr2)) {
370             len = arr2.length;
371 
372             for (i = 0; i < len; i++) {
373                 res[i] = arr1 / arr2[i];
374             }
375         } else if (Type.isArray(arr1) && Type.isArray(arr2)) {
376             len = Math.min(arr1.length, arr2.length);
377 
378             for (i = 0; i < len; i++) {
379                 res[i] = arr1[i] / arr2[i];
380             }
381         } else {
382             res = arr1 / arr2;
383         }
384 
385         return res;
386     },
387 
388     /**
389      * @function
390      * @deprecated Use {@link JXG.Math.Statistics.div} instead.
391      * @memberof JXG.Math.Statistics
392      */
393     divide: function () {
394         JXG.deprecated("Statistics.divide()", "Statistics.div()");
395         Mat.Statistics.div.apply(Mat.Statistics, arguments);
396     },
397 
398     /**
399      * Divides two (sequences of) values and returns the remainder. If two arrays are given and the lengths don't
400      * match the shortest length is taken.
401      * @param {Array|Number} arr1 Dividend
402      * @param {Array|Number} arr2 Divisor
403      * @param {Boolean} [math=false] Mathematical mod or symmetric mod? Default is symmetric, the JavaScript <tt>%</tt> operator.
404      * @returns {Array|Number}
405      * @memberof JXG.Math.Statistics
406      */
407     mod: function (arr1, arr2, math) {
408         var i,
409             len,
410             res = [],
411             mod = function (a, m) {
412                 return a % m;
413             };
414 
415         math = Type.def(math, false);
416 
417         if (math) {
418             mod = Mat.mod;
419         }
420 
421         arr1 = Type.evalSlider(arr1);
422         arr2 = Type.evalSlider(arr2);
423 
424         if (Type.isArray(arr1) && Type.isNumber(arr2)) {
425             len = arr1.length;
426 
427             for (i = 0; i < len; i++) {
428                 res[i] = mod(arr1[i], arr2);
429             }
430         } else if (Type.isNumber(arr1) && Type.isArray(arr2)) {
431             len = arr2.length;
432 
433             for (i = 0; i < len; i++) {
434                 res[i] = mod(arr1, arr2[i]);
435             }
436         } else if (Type.isArray(arr1) && Type.isArray(arr2)) {
437             len = Math.min(arr1.length, arr2.length);
438 
439             for (i = 0; i < len; i++) {
440                 res[i] = mod(arr1[i], arr2[i]);
441             }
442         } else {
443             res = mod(arr1, arr2);
444         }
445 
446         return res;
447     },
448 
449     /**
450      * Multiplies two (sequences of) values. If one value is an array and the other one is a number the number
451      * is multiplied to every element of the array. If two arrays are given and the lengths don't match the shortest
452      * length is taken.
453      * @param {Array|Number} arr1
454      * @param {Array|Number} arr2
455      * @returns {Array|Number}
456      * @memberof JXG.Math.Statistics
457      */
458     multiply: function (arr1, arr2) {
459         var i,
460             len,
461             res = [];
462 
463         arr1 = Type.evalSlider(arr1);
464         arr2 = Type.evalSlider(arr2);
465 
466         if (Type.isArray(arr1) && Type.isNumber(arr2)) {
467             len = arr1.length;
468 
469             for (i = 0; i < len; i++) {
470                 res[i] = arr1[i] * arr2;
471             }
472         } else if (Type.isNumber(arr1) && Type.isArray(arr2)) {
473             len = arr2.length;
474 
475             for (i = 0; i < len; i++) {
476                 res[i] = arr1 * arr2[i];
477             }
478         } else if (Type.isArray(arr1) && Type.isArray(arr2)) {
479             len = Math.min(arr1.length, arr2.length);
480 
481             for (i = 0; i < len; i++) {
482                 res[i] = arr1[i] * arr2[i];
483             }
484         } else {
485             res = arr1 * arr2;
486         }
487 
488         return res;
489     },
490 
491     /**
492      * Subtracts two (sequences of) values. If two arrays are given and the lengths don't match the shortest
493      * length is taken.
494      * @param {Array|Number} arr1 Minuend
495      * @param {Array|Number} arr2 Subtrahend
496      * @returns {Array|Number}
497      * @memberof JXG.Math.Statistics
498      */
499     subtract: function (arr1, arr2) {
500         var i,
501             len,
502             res = [];
503 
504         arr1 = Type.evalSlider(arr1);
505         arr2 = Type.evalSlider(arr2);
506 
507         if (Type.isArray(arr1) && Type.isNumber(arr2)) {
508             len = arr1.length;
509 
510             for (i = 0; i < len; i++) {
511                 res[i] = arr1[i] - arr2;
512             }
513         } else if (Type.isNumber(arr1) && Type.isArray(arr2)) {
514             len = arr2.length;
515 
516             for (i = 0; i < len; i++) {
517                 res[i] = arr1 - arr2[i];
518             }
519         } else if (Type.isArray(arr1) && Type.isArray(arr2)) {
520             len = Math.min(arr1.length, arr2.length);
521 
522             for (i = 0; i < len; i++) {
523                 res[i] = arr1[i] - arr2[i];
524             }
525         } else {
526             res = arr1 - arr2;
527         }
528 
529         return res;
530     },
531 
532     /**
533      * The Theil-Sen estimator can be used to determine a more robust linear regression of a set of sample
534      * points than least squares regression in {@link JXG.Math.Numerics.regressionPolynomial}.
535      *
536      * If the function should be applied to an array a of points, a the coords array can be generated with
537      * JavaScript array.map:
538      *
539      * <pre>
540      * JXG.Math.Statistics.TheilSenRegression(a.map(el => el.coords));
541      * </pre>
542      *
543      * @param {Array} coords Array of {@link JXG.Coords}.
544      * @returns {Array} A stdform array of the regression line.
545      * @memberof JXG.Math.Statistics
546      *
547      * @example
548      * var board = JXG.JSXGraph.initBoard('jxgbox', { boundingbox: [-6,6,6,-6], axis : true });
549      * var a=[];
550      * a[0]=board.create('point', [0,0]);
551      * a[1]=board.create('point', [3,0]);
552      * a[2]=board.create('point', [0,3]);
553      *
554      * board.create('line', [
555      *     () => JXG.Math.Statistics.TheilSenRegression(a.map(el => el.coords))
556      *   ],
557      *   {strokeWidth:1, strokeColor:'black'});
558      *
559      * </pre><div id="JXG0a28be85-91c5-44d3-aae6-114e81217cf0" class="jxgbox" style="width: 300px; height: 300px;"></div>
560      * <script type="text/javascript">
561      *     (function() {
562      *         var board = JXG.JSXGraph.initBoard('JXG0a28be85-91c5-44d3-aae6-114e81217cf0',
563      *             {boundingbox: [-6,6,6,-6], axis: true, showcopyright: false, shownavigation: false});
564      *     var a=[];
565      *     a[0]=board.create('point', [0,0]);
566      *     a[1]=board.create('point', [3,0]);
567      *     a[2]=board.create('point', [0,3]);
568      *
569      *     board.create('line', [
570      *         () => JXG.Math.Statistics.TheilSenRegression(a.map(el => el.coords))
571      *       ],
572      *       {strokeWidth:1, strokeColor:'black'});
573      *
574      *     })();
575      *
576      * </script><pre>
577      *
578      */
579     TheilSenRegression: function (coords) {
580         var i,
581             j,
582             slopes = [],
583             tmpslopes = [],
584             yintercepts = [];
585 
586         for (i = 0; i < coords.length; i++) {
587             tmpslopes.length = 0;
588 
589             for (j = 0; j < coords.length; j++) {
590                 if (Math.abs(coords[j].usrCoords[1] - coords[i].usrCoords[1]) > Mat.eps) {
591                     tmpslopes[j] =
592                         (coords[j].usrCoords[2] - coords[i].usrCoords[2]) /
593                         (coords[j].usrCoords[1] - coords[i].usrCoords[1]);
594                 }
595             }
596 
597             slopes[i] = this.median(tmpslopes);
598             yintercepts.push(coords[i].usrCoords[2] - slopes[i] * coords[i].usrCoords[1]);
599         }
600 
601         return [this.median(yintercepts), this.median(slopes), -1];
602     },
603 
604     /**
605      * Generate values of a standard normal random variable with the Marsaglia polar method, see
606      * {@link https://en.wikipedia.org/wiki/Marsaglia_polar_method}.
607      * See also D. E. Knuth, The art of computer programming, vol 2, p. 117.
608      *
609      * @param {Number} mean mean value of the normal distribution
610      * @param {Number} stdDev standard deviation of the normal distribution
611      * @returns {Number} value of a standard normal random variable
612      * @memberof JXG.Math.Statistics
613      */
614     generateGaussian: function (mean, stdDev) {
615         var u, v, s;
616 
617         if (this.hasSpare) {
618             this.hasSpare = false;
619             return this.spare * stdDev + mean;
620         }
621 
622         do {
623             u = Math.random() * 2 - 1;
624             v = Math.random() * 2 - 1;
625             s = u * u + v * v;
626         } while (s >= 1 || s === 0);
627 
628         s = Math.sqrt((-2.0 * Math.log(s)) / s);
629 
630         this.spare = v * s;
631         this.hasSpare = true;
632         return mean + stdDev * u * s;
633     },
634 
635     /**
636      * Generate value of a standard normal random variable with given mean and standard deviation.
637      * Alias for {@link JXG.Math.Statistics#generateGaussian}
638      *
639      * @param {Number} mean
640      * @param {Number} stdDev
641      * @returns Number
642      * @memberof JXG.Math.Statistics
643      * @see JXG.Math.Statistics#generateGaussian
644      */
645     randomNormal: function (mean, stdDev) {
646         return this.generateGaussian(mean, stdDev);
647     },
648 
649     /**
650      * Generate value of a uniform distributed random variable in the interval [a, b].
651      * @param {Number} a
652      * @param {Number} b
653      * @returns Number
654      * @memberof JXG.Math.Statistics
655      */
656     randomUniform: function (a, b) {
657         return Math.random() * (b - a) + a;
658     },
659 
660     /**
661      * Generate value of a random variable with exponential distribution, i.e.
662      * <i>f(x; lambda) = lambda * e^(-lambda x)</i> if <i>x >= 0</i> and <i>f(x; lambda) = 0</i> if <i>x < 0</i>.
663      * See {@link https://en.wikipedia.org/wiki/Exponential_distribution}.
664      * Algorithm: D.E. Knuth, TAOCP 2, p. 128.
665      *
666      * @param {Number} lambda <i>> 0</i>
667      * @returns Number
668      * @memberof JXG.Math.Statistics
669      */
670     randomExponential: function (lbda) {
671         var u;
672 
673         // Knuth, TAOCP 2, p 128
674         // See https://en.wikipedia.org/wiki/Exponential_distribution
675         if (lbda <= 0) {
676             return NaN;
677         }
678 
679         do {
680             u = Math.random();
681         } while (u === 0);
682 
683         return -Math.log(u) / lbda;
684     },
685 
686     /**
687      * Generate value of a random variable with gamma distribution of order alpha.
688      * See {@link https://en.wikipedia.org/wiki/Gamma_distribution}.
689      * Algorithm: D.E. Knuth, TAOCP 2, p. 129.
690 
691      * @param {Number} a shape, <i> > 0</i>
692      * @param {Number} [b=1] scale, <i> > 0</i>
693      * @param {Number} [t=0] threshold
694      * @returns Number
695      * @memberof JXG.Math.Statistics
696      */
697     randomGamma: function (a, b, t) {
698         var u, v, x, y,
699             p, q;
700 
701         if (a <= 0) {
702             return NaN;
703         }
704 
705         b = b || 1;
706         t = t || 0;
707 
708         if (a === 1) {
709             return b * this.randomExponential(1) + t;
710         }
711 
712         if (a < 1) {
713             // Method by Ahrens
714             // Knuth, TAOCP 2, Ex. 16, p 551
715             p = Math.E / (a + Math.E);
716 
717             do {
718                 u = Math.random();
719                 do {
720                     v = Math.random();
721                 } while (v === 0);
722                 if (u < p) {
723                     x = Math.pow(v, 1 / a);
724                     q = Math.exp(-x);
725                 } else {
726                     x = 1 - Math.log(v);
727                     q = Math.pow(x, a - 1);
728                 }
729                 u = Math.random();
730             } while (u >= q);
731             return b * x + t;
732         }
733 
734         // a > 1
735         // Knuth, TAOCP 2, p 129
736         do {
737             y = Math.tan(Math.PI * Math.random());
738             x = Math.sqrt(2 * a - 1) * y + a - 1;
739             if (x > 0) {
740                 v = Math.random();
741             } else {
742                 continue;
743             }
744         } while (x <= 0.0 || v > (1 + y * y) * Math.exp( (a - 1) * Math.log(x / (a-1)) - Math.sqrt(2 * a - 1) * y));
745 
746         return b * x + t;
747     },
748 
749     /**
750      * Generate value of a random variable with beta distribution with shape parameters alpha and beta.
751      * See {@link https://en.wikipedia.org/wiki/Beta_distribution}.
752      *
753      * @param {Number} alpha <i>> 0</i>
754      * @param {Number} beta <i>> 0</i>
755      * @returns Number
756      * @memberof JXG.Math.Statistics
757      */
758     randomBeta: function (a, b) {
759         // Knuth, TAOCP 2, p 129
760         var x1, x2, x;
761 
762         if (a <= 0 || b <= 0) {
763             return NaN;
764         }
765 
766         x1 = this.randomGamma(a);
767         x2 = this.randomGamma(b);
768         x = x1 / (x1 + x2);
769         return x;
770     },
771 
772     /**
773      * Generate value of a random variable with chi-square distribution with k degrees of freedom.
774      * See {@link https://en.wikipedia.org/wiki/Chi-squared_distribution}.
775      *
776      * @param {Number} k <i>> 0</i>
777      * @returns Number
778      * @memberof JXG.Math.Statistics
779      */
780     randomChisquare: function (nu) {
781         // Knuth, TAOCP 2, p 130
782 
783         if (nu <= 0) {
784             return NaN;
785         }
786 
787         return 2 * this.randomGamma(nu * 0.5);
788     },
789 
790     /**
791      * Generate value of a random variable with F-distribution with d<sub>1</sub> and d<sub>2</sub> degrees of freedom.
792      * See {@link https://en.wikipedia.org/wiki/F-distribution}.
793      * @param {Number} d1 <i>> 0</i>
794      * @param {Number} d2 <i>> 0</i>
795      * @returns Number
796      * @memberof JXG.Math.Statistics
797      */
798     randomF: function (nu1, nu2) {
799         // Knuth, TAOCP 2, p 130
800         var y1, y2;
801 
802         if (nu1 <= 0 || nu2 <= 0) {
803             return NaN;
804         }
805 
806         y1 = this.randomChisquare(nu1);
807         y2 = this.randomChisquare(nu2);
808 
809         return (y1 * nu2) / (y2 * nu1);
810     },
811 
812     /**
813      * Generate value of a random variable with Students-t-distribution with ν degrees of freedom.
814      * See {@link https://en.wikipedia.org/wiki/Student%27s_t-distribution}.
815      * @param {Number} nu <i>> 0</i>
816      * @returns Number
817      * @memberof JXG.Math.Statistics
818      */
819     randomT: function (nu) {
820         // Knuth, TAOCP 2, p 130
821         var y1, y2;
822 
823         if (nu <= 0) {
824             return NaN;
825         }
826 
827         y1 = this.randomNormal(0, 1);
828         y2 = this.randomChisquare(nu);
829 
830         return y1 / Math.sqrt(y2 / nu);
831     },
832 
833     /**
834      * Generate values for a random variable in binomial distribution with parameters <i>n</i> and <i>p</i>.
835      * See {@link https://en.wikipedia.org/wiki/Binomial_distribution}.
836      * It uses algorithm BG from {@link https://dl.acm.org/doi/pdf/10.1145/42372.42381}.
837      *
838      * @param {Number} n Number of trials (n >= 0)
839      * @param {Number} p Propability (0 <= p <= 1)
840      * @returns Number Integer value of a random variable in binomial distribution
841      * @memberof JXG.Math.Statistics
842      *
843      * @example
844      * console.log(JXG.Mat.Statistics.generateBinomial(100,0.1));
845      * // Possible output: 18
846      *
847      */
848     randomBinomial: function (n, p) {
849         var x, y, c,
850             a, b, N1;
851 
852         if (p < 0 || p > 1 || n < 0) {
853             return NaN;
854         }
855 
856         // Edge cases
857         if (p === 0) {
858             return 0;
859         }
860         if (p === 1) {
861             return n;
862         }
863 
864         // Now, we can assume 0 < p < 1.
865 
866         // Fast path for common cases
867         if (n === 0) {
868             return 0;
869         }
870         if (n === 1) {
871             return ((Math.random() < p) ? 1 : 0);
872         }
873 
874         // Exploit symmetry
875         if (p > 0.5) {
876             return n - this.randomBinomial(n, 1 - p);
877         }
878 
879         // General case: n > 1, p <= 0.5
880         if (n < 100) {
881             // n small:
882             // Algorithm BG (Devroye) from:
883             // https://dl.acm.org/doi/pdf/10.1145/42372.42381
884             // Time O(np) so suitable for np small only.
885             x = -1;
886             y = 0;
887 
888             c = Math.log(1 - p);
889             if (c === 0) {
890                 return 0;
891             }
892 
893             do {
894                 x += 1;
895                 y += Math.floor(Math.log(Math.random()) / c) + 1;
896             } while (y < n);
897         } else {
898             // n large:
899             // Knuth, TAOCP 2, p 131
900             a = 1 + Math.floor(n * 0.5);
901             b = n - a + 1;
902             x = this.randomBeta(a, b);
903             if (x >= p) {
904                 N1 = this.randomBinomial(a - 1, p / x);
905                 x = N1;
906             } else {
907                 N1 = this.randomBinomial(b - 1, (p - x) / (1 - x));
908                 x = a + N1;
909             }
910         }
911         return x;
912     },
913 
914     /**
915      * Generate values for a random variable in geometric distribution with propability <i>p</i>.
916      * See {@link https://en.wikipedia.org/wiki/Geometric_distribution}.
917      *
918      * @param {Number} p (0 <= p <= 1)
919      * @returns Number
920      * @memberof JXG.Math.Statistics
921      */
922     randomGeometric: function(p) {
923         var u;
924 
925         if (p < 0 || p > 1) {
926             return NaN;
927         }
928         // Knuth, TAOCP 2, p 131
929         u = Math.random();
930 
931         return Math.ceil(Math.log(u) / Math.log(1 - p));
932     },
933 
934     /**
935      * Generate values for a random variable in Poisson distribution with mean <i>mu</i>.
936      * See {@link https://en.wikipedia.org/wiki/Poisson_distribution}.
937      *
938      * @param {Number} mu (0 < mu)
939      * @returns Number
940      * @memberof JXG.Math.Statistics
941      */
942     randomPoisson: function (mu) {
943         var e = Math.exp(-mu),
944             N,
945             m = 0,
946             u = 1,
947             x,
948             alpha = 7 / 8;
949 
950         if (mu <= 0) {
951             return NaN;
952         }
953 
954         // Knuth, TAOCP 2, p 132
955         if (mu < 10) {
956             do {
957                 u *= Math.random();
958                 m += 1;
959             } while (u > e);
960             N = m - 1;
961         } else {
962             m = Math.floor(alpha * mu);
963             x = this.randomGamma(m);
964             if (x < mu) {
965                 N = m + this.randomPoisson(mu - x);
966             } else {
967                 N = this.randomBinomial(m - 1, mu / x);
968             }
969         }
970         return N;
971     },
972 
973     /**
974      * Generate values for a random variable in hypergeometric distribution.
975      * Samples are drawn from a hypergeometric distribution with specified parameters, <i>good</i> (ways to make a good selection),
976      * <i>bad</i> (ways to make a bad selection), and <i>samples</i> (number of items sampled, which is less than or equal to <i>good + bad</i>).
977      * <p>
978      * Naive implementation with runtime <i>O(samples)</i>.
979      *
980      * @param {Number} good ways to make a good selection
981      * @param {Number} bad ways to make a bad selection
982      * @param {Number} samples number of items sampled
983      * @returns
984      * @memberof JXG.Math.Statistics
985      */
986     randomHypergeometric: function (good, bad, k) {
987         var i, u,
988             x = 0,
989             // kk,
990             // n = good + bad,
991             d1 = good + bad - k,
992             d2 = Math.min(good, bad),
993             y = d2;
994 
995         if (good < 1 || bad < 1 || k > good + bad) {
996             return NaN;
997         }
998 
999         // Naive method
1000         // kk = Math.min(k, n - k);
1001         // for (i = 0; i < k; i ++) {
1002         //     u = Math.random();
1003         //     if (n * u <= good) {
1004         //         x += 1;
1005         //         if (x === good) {
1006         //             return x;
1007         //         }
1008         //         good -= 1;
1009         //     }
1010         //     n -= 1;
1011         // }
1012         // return x;
1013 
1014         // Implementation from
1015         // Monte Carlo by George S. Fishman
1016         // https://link.springer.com/book/10.1007/978-1-4757-2553-7
1017         // page 218
1018         //
1019         i = k;
1020         while (y * i > 0) {
1021             u = Math.random();
1022             y -= Math.floor(u + y / (d1 + i));
1023             i -= 1;
1024         }
1025         x = d2 - y;
1026         if (good <= bad) {
1027             return x;
1028         } else {
1029             return k - x;
1030         }
1031     },
1032 
1033     /**
1034      * Compute the histogram of a dataset.
1035      * Optional parameters can be supplied through a JavaScript object
1036      * with the following default values:
1037      * <pre>
1038      * {
1039      *   bins: 10,          // Number of bins
1040      *   range: false,      // false or array. The lower and upper range of the bins.
1041      *                      // If not provided, range is simply [min(x), max(x)].
1042      *                      // Values outside the range are ignored.
1043      *   density: false,    // If true, normalize the counts by dividing by sum(counts)
1044      *   cumulative: false
1045      * }
1046      * </pre>
1047      * The function returns an array containing two arrays. The first array is of length bins+1
1048      * containing the start values of the bins. The last entry contains the end values of the last bin.
1049      * <p>
1050      * The second array contains the counts of each bin.
1051      * @param {Array} x
1052      * @param {Object} opt Optional parameters
1053      * @returns Array [bin, counts] Array bins contains start values of bins, array counts contains
1054      * the number of entries of x which are contained in each bin.
1055      * @memberof JXG.Math.Statistics
1056      *
1057      * @example
1058      *     var curve = board.create('curve', [[], []]);
1059      *     curve.updateDataArray = function () {
1060      *       var i, res, x = [];
1061      *
1062      *       for (i = 0; i < 5000; i++) {
1063      *         x.push(JXG.Math.Statistics.randomGamma(2));
1064      *       }
1065      *       res = JXG.Math.Statistics.histogram(x, { bins: 50, density: true, cumulative: false, range: [-5, 5] });
1066      *       this.dataX = res[1];
1067      *       this.dataY = res[0];
1068      *     };
1069      *     board.update();
1070      *
1071      * </pre><div id="JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302" class="jxgbox" style="width: 300px; height: 300px;"></div>
1072      * <script type="text/javascript">
1073      *     (function() {
1074      *         var board = JXG.JSXGraph.initBoard('JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302',
1075      *             {boundingbox: [-1, 3, 6, -1], axis: true, showcopyright: false, shownavigation: false});
1076      *         var curve = board.create('curve', [[], []]);
1077      *         curve.updateDataArray = function () {
1078      *           var i, res, x = [];
1079      *
1080      *           for (i = 0; i < 5000; i++) {
1081      *             x.push(JXG.Math.Statistics.randomGamma(2));
1082      *           }
1083      *           res = JXG.Math.Statistics.histogram(x, { bins: 50, density: true, cumulative: false, range: [-5, 5] });
1084      *           this.dataX = res[1];
1085      *           this.dataY = res[0];
1086      *         };
1087      *         board.update();
1088      *     })();
1089      *
1090      * </script><pre>
1091      *
1092      */
1093     histogram: function (x, opt) {
1094         var i, le, k,
1095             mi, ma, num_bins, delta,
1096             range,
1097             s,
1098             counts = [],
1099             bins = [];
1100 
1101         // Evaluate number of bins
1102         num_bins = opt.bins || 10;
1103 
1104         // Evaluate range
1105         range = opt.range || false;
1106         if (range === false) {
1107             mi = Math.min.apply(null, x);
1108             ma = Math.max.apply(null, x);
1109         } else {
1110             mi = range[0];
1111             ma = range[1];
1112         }
1113 
1114         // Set uniform delta
1115         if (num_bins > 0) {
1116             delta = (ma - mi) / (num_bins - 1);
1117         } else {
1118             delta = 0;
1119         }
1120 
1121         // Set the bins and init the counts array
1122         for (i = 0; i < num_bins; i++) {
1123             counts.push(0);
1124             bins.push(mi + i * delta);
1125         }
1126         // bins.push(ma);
1127 
1128         // Determine the counts
1129         le = x.length;
1130         for (i = 0; i < le; i++) {
1131             k = Math.floor((x[i] - mi) / delta);
1132             if (k >= 0 && k < num_bins) {
1133                 counts[k] += 1;
1134             }
1135         }
1136 
1137         // Normalize if density===true
1138         if (opt.density) {
1139             s = JXG.Math.Statistics.sum(counts);
1140             for (i = 0; i < num_bins; i++) {
1141                 // counts[i] /= (s * delta);
1142                 counts[i] /= s;
1143             }
1144         }
1145 
1146         // Cumulative counts
1147         if (opt.cumulative) {
1148             for (i = 1; i < num_bins; i++) {
1149                 counts[i] += counts[i - 1];
1150             }
1151         }
1152 
1153         return [counts, bins];
1154     }
1155 };
1156 
1157 export default Mat.Statistics;
1158