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 ( <i>0 < P ≤ 100</i> ) of a list of <i>N</i> ordered values (sorted from least to greatest)
129      * is the smallest value in the list such that no more than <i>P</i> percent of the data is strictly less
130      * than the value and at least <i>P</i> percent of the data is less than or equal to that value.
131      * See <a href="https://en.wikipedia.org/wiki/Percentile">https://en.wikipedia.org/wiki/Percentile</a>.
132      *
133      * Here, the <i>linear interpolation between closest ranks</i> method is used.
134      * @param {Array} arr The set of values, need not be ordered.
135      * @param {Number|Array} percentile One or several percentiles
136      * @returns {Number|Array} Depending if a number or an array is the input for percentile, a number or an array containing the percentiles
137      * is returned.
138      */
139     percentile: function (arr, percentile) {
140         var tmp,
141             len,
142             i,
143             p,
144             res = [],
145             per;
146 
147         if (arr.length > 0) {
148             if (ArrayBuffer.isView(arr)) {
149                 tmp = new Float64Array(arr);
150                 tmp.sort();
151             } else {
152                 tmp = arr.slice(0);
153                 tmp.sort(function (a, b) {
154                     return a - b;
155                 });
156             }
157             len = tmp.length;
158 
159             if (Type.isArray(percentile)) {
160                 p = percentile;
161             } else {
162                 p = [percentile];
163             }
164 
165             for (i = 0; i < p.length; i++) {
166                 per = len * p[i] * 0.01;
167                 if (parseInt(per, 10) === per) {
168                     res.push((tmp[per - 1] + tmp[per]) * 0.5);
169                 } else {
170                     res.push(tmp[parseInt(per, 10)]);
171                 }
172             }
173 
174             if (Type.isArray(percentile)) {
175                 return res;
176             } else {
177                 return res[0];
178             }
179         }
180 
181         return 0.0;
182     },
183 
184     /**
185      * Bias-corrected sample variance. A variance is a measure of how far a
186      * set of numbers are spread out from each other.
187      * @param {Array} arr
188      * @returns {Number}
189      * @memberof JXG.Math.Statistics
190      */
191     variance: function (arr) {
192         var m,
193             res,
194             i,
195             len = arr.length;
196 
197         if (len > 1) {
198             m = this.mean(arr);
199             res = 0;
200             for (i = 0; i < len; i++) {
201                 res += (arr[i] - m) * (arr[i] - m);
202             }
203             return res / (arr.length - 1);
204         }
205 
206         return 0.0;
207     },
208 
209     /**
210      * Determines the <strong>s</strong>tandard <strong>d</strong>eviation which shows how much
211      * variation there is from the average value of a set of numbers.
212      * @param {Array} arr
213      * @returns {Number}
214      * @memberof JXG.Math.Statistics
215      */
216     sd: function (arr) {
217         return Math.sqrt(this.variance(arr));
218     },
219 
220     /**
221      * Weighted mean value is basically the same as {@link JXG.Math.Statistics.mean} but here the values
222      * are weighted, i.e. multiplied with another value called <em>weight</em>. The weight values are given
223      * as a second array with the same length as the value array..
224      * @throws {Error} If the dimensions of the arrays don't match.
225      * @param {Array} arr Set of alues.
226      * @param {Array} w Weight values.
227      * @returns {Number}
228      * @memberof JXG.Math.Statistics
229      */
230     weightedMean: function (arr, w) {
231         if (arr.length !== w.length) {
232             throw new Error(
233                 "JSXGraph error (Math.Statistics.weightedMean): Array dimension mismatch."
234             );
235         }
236 
237         if (arr.length > 0) {
238             return this.mean(this.multiply(arr, w));
239         }
240 
241         return 0.0;
242     },
243 
244     /**
245      * Extracts the maximum value from the array.
246      * @param {Array} arr
247      * @returns {Number} The highest number from the array. It returns <tt>NaN</tt> if not every element could be
248      * interpreted as a number and <tt>-Infinity</tt> if an empty array is given or no element could be interpreted
249      * as a number.
250      * @memberof JXG.Math.Statistics
251      */
252     max: function (arr) {
253         return Math.max.apply(this, arr);
254     },
255 
256     /**
257      * Extracts the minimum value from the array.
258      * @param {Array} arr
259      * @returns {Number} The lowest number from the array. It returns <tt>NaN</tt> if not every element could be
260      * interpreted as a number and <tt>Infinity</tt> if an empty array is given or no element could be interpreted
261      * as a number.
262      * @memberof JXG.Math.Statistics
263      */
264     min: function (arr) {
265         return Math.min.apply(this, arr);
266     },
267 
268     /**
269      * Determines the lowest and the highest value from the given array.
270      * @param {Array} arr
271      * @returns {Array} The minimum value as the first and the maximum value as the second value.
272      * @memberof JXG.Math.Statistics
273      */
274     range: function (arr) {
275         return [this.min(arr), this.max(arr)];
276     },
277 
278     /**
279      * Determines the absolute value of every given value.
280      * @param {Array|Number} arr
281      * @returns {Array|Number}
282      * @memberof JXG.Math.Statistics
283      */
284     abs: function (arr) {
285         var i, len, res;
286 
287         if (Type.isArray(arr)) {
288             if (arr.map) {
289                 res = arr.map(Math.abs);
290             } else {
291                 len = arr.length;
292                 res = [];
293 
294                 for (i = 0; i < len; i++) {
295                     res[i] = Math.abs(arr[i]);
296                 }
297             }
298         } else if (ArrayBuffer.isView(arr)) {
299             res = arr.map(Math.abs);
300         } else {
301             res = Math.abs(arr);
302         }
303         return res;
304     },
305 
306     /**
307      * Adds up two (sequences of) values. If one value is an array and the other one is a number the number
308      * is added to every element of the array. If two arrays are given and the lengths don't match the shortest
309      * length is taken.
310      * @param {Array|Number} arr1
311      * @param {Array|Number} arr2
312      * @returns {Array|Number}
313      * @memberof JXG.Math.Statistics
314      */
315     add: function (arr1, arr2) {
316         var i,
317             len,
318             res = [];
319 
320         arr1 = Type.evalSlider(arr1);
321         arr2 = Type.evalSlider(arr2);
322 
323         if (Type.isArray(arr1) && Type.isNumber(arr2)) {
324             len = arr1.length;
325 
326             for (i = 0; i < len; i++) {
327                 res[i] = arr1[i] + arr2;
328             }
329         } else if (Type.isNumber(arr1) && Type.isArray(arr2)) {
330             len = arr2.length;
331 
332             for (i = 0; i < len; i++) {
333                 res[i] = arr1 + arr2[i];
334             }
335         } else if (Type.isArray(arr1) && Type.isArray(arr2)) {
336             len = Math.min(arr1.length, arr2.length);
337 
338             for (i = 0; i < len; i++) {
339                 res[i] = arr1[i] + arr2[i];
340             }
341         } else {
342             res = arr1 + arr2;
343         }
344 
345         return res;
346     },
347 
348     /**
349      * Divides two (sequences of) values. If two arrays are given and the lengths don't match the shortest length
350      * is taken.
351      * @param {Array|Number} arr1 Dividend
352      * @param {Array|Number} arr2 Divisor
353      * @returns {Array|Number}
354      * @memberof JXG.Math.Statistics
355      */
356     div: function (arr1, arr2) {
357         var i,
358             len,
359             res = [];
360 
361         arr1 = Type.evalSlider(arr1);
362         arr2 = Type.evalSlider(arr2);
363 
364         if (Type.isArray(arr1) && Type.isNumber(arr2)) {
365             len = arr1.length;
366 
367             for (i = 0; i < len; i++) {
368                 res[i] = arr1[i] / arr2;
369             }
370         } else if (Type.isNumber(arr1) && Type.isArray(arr2)) {
371             len = arr2.length;
372 
373             for (i = 0; i < len; i++) {
374                 res[i] = arr1 / arr2[i];
375             }
376         } else if (Type.isArray(arr1) && Type.isArray(arr2)) {
377             len = Math.min(arr1.length, arr2.length);
378 
379             for (i = 0; i < len; i++) {
380                 res[i] = arr1[i] / arr2[i];
381             }
382         } else {
383             res = arr1 / arr2;
384         }
385 
386         return res;
387     },
388 
389     /**
390      * @function
391      * @deprecated Use {@link JXG.Math.Statistics.div} instead.
392      * @memberof JXG.Math.Statistics
393      */
394     divide: function () {
395         JXG.deprecated("Statistics.divide()", "Statistics.div()");
396         Mat.Statistics.div.apply(Mat.Statistics, arguments);
397     },
398 
399     /**
400      * Divides two (sequences of) values and returns the remainder. If two arrays are given and the lengths don't
401      * match the shortest length is taken.
402      * @param {Array|Number} arr1 Dividend
403      * @param {Array|Number} arr2 Divisor
404      * @param {Boolean} [math=false] Mathematical mod or symmetric mod? Default is symmetric, the JavaScript <tt>%</tt> operator.
405      * @returns {Array|Number}
406      * @memberof JXG.Math.Statistics
407      */
408     mod: function (arr1, arr2, math) {
409         var i,
410             len,
411             res = [],
412             mod = function (a, m) {
413                 return a % m;
414             };
415 
416         math = Type.def(math, false);
417 
418         if (math) {
419             mod = Mat.mod;
420         }
421 
422         arr1 = Type.evalSlider(arr1);
423         arr2 = Type.evalSlider(arr2);
424 
425         if (Type.isArray(arr1) && Type.isNumber(arr2)) {
426             len = arr1.length;
427 
428             for (i = 0; i < len; i++) {
429                 res[i] = mod(arr1[i], arr2);
430             }
431         } else if (Type.isNumber(arr1) && Type.isArray(arr2)) {
432             len = arr2.length;
433 
434             for (i = 0; i < len; i++) {
435                 res[i] = mod(arr1, arr2[i]);
436             }
437         } else if (Type.isArray(arr1) && Type.isArray(arr2)) {
438             len = Math.min(arr1.length, arr2.length);
439 
440             for (i = 0; i < len; i++) {
441                 res[i] = mod(arr1[i], arr2[i]);
442             }
443         } else {
444             res = mod(arr1, arr2);
445         }
446 
447         return res;
448     },
449 
450     /**
451      * Multiplies two (sequences of) values. If one value is an array and the other one is a number the number
452      * is multiplied to every element of the array. If two arrays are given and the lengths don't match the shortest
453      * length is taken.
454      * @param {Array|Number} arr1
455      * @param {Array|Number} arr2
456      * @returns {Array|Number}
457      * @memberof JXG.Math.Statistics
458      */
459     multiply: function (arr1, arr2) {
460         var i,
461             len,
462             res = [];
463 
464         arr1 = Type.evalSlider(arr1);
465         arr2 = Type.evalSlider(arr2);
466 
467         if (Type.isArray(arr1) && Type.isNumber(arr2)) {
468             len = arr1.length;
469 
470             for (i = 0; i < len; i++) {
471                 res[i] = arr1[i] * arr2;
472             }
473         } else if (Type.isNumber(arr1) && Type.isArray(arr2)) {
474             len = arr2.length;
475 
476             for (i = 0; i < len; i++) {
477                 res[i] = arr1 * arr2[i];
478             }
479         } else if (Type.isArray(arr1) && Type.isArray(arr2)) {
480             len = Math.min(arr1.length, arr2.length);
481 
482             for (i = 0; i < len; i++) {
483                 res[i] = arr1[i] * arr2[i];
484             }
485         } else {
486             res = arr1 * arr2;
487         }
488 
489         return res;
490     },
491 
492     /**
493      * Subtracts two (sequences of) values. If two arrays are given and the lengths don't match the shortest
494      * length is taken.
495      * @param {Array|Number} arr1 Minuend
496      * @param {Array|Number} arr2 Subtrahend
497      * @returns {Array|Number}
498      * @memberof JXG.Math.Statistics
499      */
500     subtract: function (arr1, arr2) {
501         var i,
502             len,
503             res = [];
504 
505         arr1 = Type.evalSlider(arr1);
506         arr2 = Type.evalSlider(arr2);
507 
508         if (Type.isArray(arr1) && Type.isNumber(arr2)) {
509             len = arr1.length;
510 
511             for (i = 0; i < len; i++) {
512                 res[i] = arr1[i] - arr2;
513             }
514         } else if (Type.isNumber(arr1) && Type.isArray(arr2)) {
515             len = arr2.length;
516 
517             for (i = 0; i < len; i++) {
518                 res[i] = arr1 - arr2[i];
519             }
520         } else if (Type.isArray(arr1) && Type.isArray(arr2)) {
521             len = Math.min(arr1.length, arr2.length);
522 
523             for (i = 0; i < len; i++) {
524                 res[i] = arr1[i] - arr2[i];
525             }
526         } else {
527             res = arr1 - arr2;
528         }
529 
530         return res;
531     },
532 
533     /**
534      * The Theil-Sen estimator can be used to determine a more robust linear regression of a set of sample
535      * points than least squares regression in {@link JXG.Math.Numerics.regressionPolynomial}.
536      *
537      * If the function should be applied to an array a of points, a the coords array can be generated with
538      * JavaScript array.map:
539      *
540      * <pre>
541      * JXG.Math.Statistics.TheilSenRegression(a.map(el => el.coords));
542      * </pre>
543      *
544      * @param {Array} coords Array of {@link JXG.Coords}.
545      * @returns {Array} A stdform array of the regression line.
546      * @memberof JXG.Math.Statistics
547      *
548      * @example
549      * var board = JXG.JSXGraph.initBoard('jxgbox', { boundingbox: [-6,6,6,-6], axis : true });
550      * var a=[];
551      * a[0]=board.create('point', [0,0]);
552      * a[1]=board.create('point', [3,0]);
553      * a[2]=board.create('point', [0,3]);
554      *
555      * board.create('line', [
556      *     () => JXG.Math.Statistics.TheilSenRegression(a.map(el => el.coords))
557      *   ],
558      *   {strokeWidth:1, strokeColor:'black'});
559      *
560      * </pre><div id="JXG0a28be85-91c5-44d3-aae6-114e81217cf0" class="jxgbox" style="width: 300px; height: 300px;"></div>
561      * <script type="text/javascript">
562      *     (function() {
563      *         var board = JXG.JSXGraph.initBoard('JXG0a28be85-91c5-44d3-aae6-114e81217cf0',
564      *             {boundingbox: [-6,6,6,-6], axis: true, showcopyright: false, shownavigation: false});
565      *     var a=[];
566      *     a[0]=board.create('point', [0,0]);
567      *     a[1]=board.create('point', [3,0]);
568      *     a[2]=board.create('point', [0,3]);
569      *
570      *     board.create('line', [
571      *         () => JXG.Math.Statistics.TheilSenRegression(a.map(el => el.coords))
572      *       ],
573      *       {strokeWidth:1, strokeColor:'black'});
574      *
575      *     })();
576      *
577      * </script><pre>
578      *
579      */
580     TheilSenRegression: function (coords) {
581         var i,
582             j,
583             slopes = [],
584             tmpslopes = [],
585             yintercepts = [];
586 
587         for (i = 0; i < coords.length; i++) {
588             tmpslopes.length = 0;
589 
590             for (j = 0; j < coords.length; j++) {
591                 if (Math.abs(coords[j].usrCoords[1] - coords[i].usrCoords[1]) > Mat.eps) {
592                     tmpslopes[j] =
593                         (coords[j].usrCoords[2] - coords[i].usrCoords[2]) /
594                         (coords[j].usrCoords[1] - coords[i].usrCoords[1]);
595                 }
596             }
597 
598             slopes[i] = this.median(tmpslopes);
599             yintercepts.push(coords[i].usrCoords[2] - slopes[i] * coords[i].usrCoords[1]);
600         }
601 
602         return [this.median(yintercepts), this.median(slopes), -1];
603     },
604 
605     /**
606      * Generate values of a standard normal random variable with the Marsaglia polar method, see
607      * <a href="https://en.wikipedia.org/wiki/Marsaglia_polar_method">https://en.wikipedia.org/wiki/Marsaglia_polar_method</a>.
608      * See also D. E. Knuth, The art of computer programming, vol 2, p. 117.
609      *
610      * @param {Number} mean mean value of the normal distribution
611      * @param {Number} stdDev standard deviation of the normal distribution
612      * @returns {Number} value of a standard normal random variable
613      * @memberof JXG.Math.Statistics
614      */
615     generateGaussian: function (mean, stdDev) {
616         var u, v, s;
617 
618         if (this.hasSpare) {
619             this.hasSpare = false;
620             return this.spare * stdDev + mean;
621         }
622 
623         do {
624             u = Math.random() * 2 - 1;
625             v = Math.random() * 2 - 1;
626             s = u * u + v * v;
627         } while (s >= 1 || s === 0);
628 
629         s = Math.sqrt((-2.0 * Math.log(s)) / s);
630 
631         this.spare = v * s;
632         this.hasSpare = true;
633         return mean + stdDev * u * s;
634     },
635 
636     /**
637      * Generate value of a standard normal random variable with given mean and standard deviation.
638      * Alias for {@link JXG.Math.Statistics#generateGaussian}
639      *
640      * @param {Number} mean
641      * @param {Number} stdDev
642      * @returns Number
643      * @memberof JXG.Math.Statistics
644      * @see JXG.Math.Statistics.generateGaussian
645      * @example
646      *  let board = JXG.JSXGraph.initBoard('JXGbox',
647      *       { boundingbox: [-5, 1.5, 5, -.03], axis: true});
648      *
649      *   let runs = [
650      *       [0, 0.2, 'blue'],
651      *       [0, 1.0, 'red'],
652      *       [0, 5.0, 'orange'],
653      *       [-2,0.5, 'green'],
654      *   ]
655      *
656      *   let labelY = 1.2
657      *   runs.forEach((run,i) => {
658      *       board.create('segment',[[1.0,labelY-(i/20)],[2.0,labelY-(i/20)]],{strokeColor:run[2]})
659      *       board.create('text',[2.5,labelY-(i/20),`μ=${run[0]}, σ<sup>2</sup>=${run[1]}`])
660      *
661      *       let x = Array(50000).fill(0).map(() => JXG.Math.Statistics.randomNormal(run[0],Math.sqrt(run[1])))  // sqrt so Std Dev, not Variance
662      *       let res = JXG.Math.Statistics.histogram(x, { bins: 40, density: true, cumulative: false, range: false });
663      *       board.create('curve', [res[1], res[0]], { strokeColor: run[2], strokeWidth:2});
664      *   })
665      *
666      * </pre><div id="JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302-4" class="jxgbox" style="width: 300px; height: 300px;"></div>
667      * <script type="text/javascript">
668      * {
669      *  let board = JXG.JSXGraph.initBoard('JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302-4',
670      *       { boundingbox: [-5, 1.5, 5, -.03], axis: true});
671      *
672      *   let runs = [
673      *       [0, 0.2, 'blue'],
674      *       [0, 1.0, 'red'],
675      *       [0, 5.0, 'orange'],
676      *       [-2,0.5, 'green'],
677      *   ]
678      *
679      *   let labelY = 1.2
680      *   runs.forEach((run,i) => {
681      *       board.create('segment',[[1.0,labelY-(i/20)],[2.0,labelY-(i/20)]],{strokeColor:run[2]})
682      *       board.create('text',[2.5,labelY-(i/20),`μ=${run[0]}, σ<sup>2</sup>=${run[1]}`])
683      *
684      *       let x = Array(50000).fill(0).map(() => JXG.Math.Statistics.randomNormal(run[0],Math.sqrt(run[1])))  // sqrt so Std Dev, not Variance
685      *       let res = JXG.Math.Statistics.histogram(x, { bins: 40, density: true, cumulative: false, range: false });
686      *       board.create('curve', [res[1], res[0]], { strokeColor: run[2], strokeWidth:2});
687      *   })
688      * }
689      * </script><pre>
690 
691      */
692     randomNormal: function (mean, stdDev) {
693         return this.generateGaussian(mean, stdDev);
694     },
695 
696     /**
697      * Generate value of a uniform distributed random variable in the interval [a, b].
698      * @param {Number} a
699      * @param {Number} b
700      * @returns Number
701      * @memberof JXG.Math.Statistics
702      */
703     randomUniform: function (a, b) {
704         return Math.random() * (b - a) + a;
705     },
706 
707     /**
708      * Generate value of a random variable with exponential distribution, i.e.
709      * <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>.
710      * See <a href="https://en.wikipedia.org/wiki/Exponential_distribution">https://en.wikipedia.org/wiki/Exponential_distribution</a>.
711      * Algorithm: D.E. Knuth, TAOCP 2, p. 128.
712      *
713      * @param {Number} lambda <i>> 0</i>
714      * @returns Number
715      * @memberof JXG.Math.Statistics
716      * @example
717      *  let board = JXG.JSXGraph.initBoard('JXGbox',
718      *       { boundingbox: [-.5, 1.5, 5, -.1], axis: true});
719      *
720      *   let runs = [
721      *       [0.5, 'red'],
722      *       [1.0, 'green'],
723      *       [1.5, 'blue'],
724      *   ]
725      *
726      *   let labelY = 1
727      *   runs.forEach((run,i) => {
728      *       board.create('segment',[[1.8,labelY-(i/20)],[2.3,labelY-(i/20)]],{strokeColor:run[1]})
729      *       board.create('text',[2.5,labelY-(i/20),`λ=${run[0]}`])
730      *
731      *       let x = Array(50000).fill(0).map(() => JXG.Math.Statistics.randomExponential(run[0]))
732      *       let res = JXG.Math.Statistics.histogram(x, { bins: 40, density: true, cumulative: false, range: false });
733      *       board.create('curve', [res[1], res[0]], { strokeColor: run[1], strokeWidth:2});
734      *   })
735      *
736      * </pre><div id="JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302-5" class="jxgbox" style="width: 300px; height: 300px;"></div>
737      * <script type="text/javascript">
738      * {
739      *  let board = JXG.JSXGraph.initBoard('JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302-5',
740      *       { boundingbox: [-.5, 1.5, 5, -.1], axis: true});
741      *
742      *   let runs = [
743      *       [0.5, 'red'],
744      *       [1.0, 'green'],
745      *       [1.5, 'blue'],
746      *   ]
747      *
748      *   let labelY = 1
749      *   runs.forEach((run,i) => {
750      *       board.create('segment',[[1.8,labelY-(i/20)],[2.3,labelY-(i/20)]],{strokeColor:run[1]})
751      *       board.create('text',[2.5,labelY-(i/20),`λ=${run[0]}`])
752      *
753      *       let x = Array(50000).fill(0).map(() => JXG.Math.Statistics.randomExponential(run[0]))
754      *       let res = JXG.Math.Statistics.histogram(x, { bins: 40, density: true, cumulative: false, range: false });
755      *       board.create('curve', [res[1], res[0]], { strokeColor: run[1], strokeWidth:2});
756      *   })
757      * }
758      * </script><pre>
759 
760     */
761     randomExponential: function (lbda) {
762         var u;
763 
764         // Knuth, TAOCP 2, p 128
765         // See https://en.wikipedia.org/wiki/Exponential_distribution
766         if (lbda <= 0) {
767             return NaN;
768         }
769 
770         do {
771             u = Math.random();
772         } while (u === 0);
773 
774         return -Math.log(u) / lbda;
775     },
776 
777     /**
778      * Generate value of a random variable with gamma distribution of order alpha.
779      * See <a href="https://en.wikipedia.org/wiki/Gamma_distribution">https://en.wikipedia.org/wiki/Gamma_distribution</a>.
780      * Algorithm: D.E. Knuth, TAOCP 2, p. 129.
781 
782      * @param {Number} a shape, <i> > 0</i>
783      * @param {Number} [b=1] scale, <i> > 0</i>
784      * @param {Number} [t=0] threshold
785      * @returns Number
786      * @memberof JXG.Math.Statistics
787      * @example
788      *  let board = JXG.JSXGraph.initBoard('jxgbox',
789      *       { boundingbox: [-1.7, .5, 20, -.03], axis: true});
790      *
791      *   let runs = [
792      *       [0.5, 1.0, 'brown'],
793      *       [1.0, 2.0, 'red'],
794      *       [2.0, 2.0, 'orange'],
795      *       [3.0, 2.0, 'yellow'],
796      *       [5.0, 1.0, 'green'],
797      *       [9.0, 0.5, 'black'],
798      *       [7.5, 1.0, 'purple'],
799      *   ]
800      *
801      *   let labelY = .4
802      *   runs.forEach((run,i) => {
803      *       board.create('segment',[[7,labelY-(i/50)],[9,labelY-(i/50)]],{strokeColor:run[2]})
804      *       board.create('text',[10,labelY-(i/50),`k=${run[0]}, θ=${run[1]}`])
805      *
806      *       // density
807      *       let x = Array(50000).fill(0).map(() => JXG.Math.Statistics.randomGamma(run[0],run[1]))
808      *       let res = JXG.Math.Statistics.histogram(x, { bins: 50, density: true, cumulative: false, range: [0, 20] });
809      *       board.create('curve', [res[1], res[0]], { strokeColor: run[2]});
810      *
811      *   })
812      *
813      *
814      * </pre>
815      * <div id="JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302-6" class="jxgbox" style="width: 300px; height: 300px;"></div>
816      * <script type="text/javascript">
817      * {
818      *  let board = JXG.JSXGraph.initBoard('JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302-6',
819      *       { boundingbox: [-1.7, .5, 20, -.03], axis: true});
820      *
821      *   let runs = [
822      *       [0.5, 1.0, 'brown'],
823      *       [1.0, 2.0, 'red'],
824      *       [2.0, 2.0, 'orange'],
825      *       [3.0, 2.0, 'yellow'],
826      *       [5.0, 1.0, 'green'],
827      *       [9.0, 0.5, 'black'],
828      *       [7.5, 1.0, 'purple'],
829      *   ]
830      *
831      *   let labelY = .4
832      *   runs.forEach((run,i) => {
833      *       board.create('segment',[[7,labelY-(i/50)],[9,labelY-(i/50)]],{strokeColor:run[2]})
834      *       board.create('text',[10,labelY-(i/50),`k=${run[0]}, θ=${run[1]}`])
835      *
836      *       let x = Array(50000).fill(0).map(() => JXG.Math.Statistics.randomGamma(run[0],run[1]))
837      *       let res = JXG.Math.Statistics.histogram(x, { bins: 50, density: true, cumulative: false, range: [0, 20] });
838      *       board.create('curve', [res[1], res[0]], { strokeColor: run[2]});
839      *   })
840      * }
841      * </script><pre>
842      *
843      */
844     randomGamma: function (a, b, t) {
845         var u, v, x, y,
846             p, q;
847 
848         if (a <= 0) {
849             return NaN;
850         }
851 
852         b = b || 1;
853         t = t || 0;
854 
855         if (a === 1) {
856             return b * this.randomExponential(1) + t;
857         }
858 
859         if (a < 1) {
860             // Method by Ahrens
861             // Knuth, TAOCP 2, Ex. 16, p 551
862             p = Math.E / (a + Math.E);
863 
864             do {
865                 u = Math.random();
866                 do {
867                     v = Math.random();
868                 } while (v === 0);
869                 if (u < p) {
870                     x = Math.pow(v, 1 / a);
871                     q = Math.exp(-x);
872                 } else {
873                     x = 1 - Math.log(v);
874                     q = Math.pow(x, a - 1);
875                 }
876                 u = Math.random();
877             } while (u >= q);
878             return b * x + t;
879         }
880 
881         // a > 1
882         // Knuth, TAOCP 2, p 129
883         do {
884             y = Math.tan(Math.PI * Math.random());
885             x = Math.sqrt(2 * a - 1) * y + a - 1;
886             if (x > 0) {
887                 v = Math.random();
888             } else {
889                 continue;
890             }
891         } while (x <= 0.0 || v > (1 + y * y) * Math.exp((a - 1) * Math.log(x / (a - 1)) - Math.sqrt(2 * a - 1) * y));
892 
893         return b * x + t;
894     },
895 
896     /**
897      * Generate value of a random variable with beta distribution with shape parameters alpha and beta.
898      * See <a href="https://en.wikipedia.org/wiki/Beta_distribution">https://en.wikipedia.org/wiki/Beta_distribution</a>.
899      *
900      * @param {Number} alpha <i>> 0</i>
901      * @param {Number} beta <i>> 0</i>
902      * @returns Number
903      * @memberof JXG.Math.Statistics
904      */
905     randomBeta: function (a, b) {
906         // Knuth, TAOCP 2, p 129
907         var x1, x2, x;
908 
909         if (a <= 0 || b <= 0) {
910             return NaN;
911         }
912 
913         x1 = this.randomGamma(a);
914         x2 = this.randomGamma(b);
915         x = x1 / (x1 + x2);
916         return x;
917     },
918 
919     /**
920      * Generate value of a random variable with chi-square distribution with k degrees of freedom.
921      * See <a href="https://en.wikipedia.org/wiki/Chi-squared_distribution">https://en.wikipedia.org/wiki/Chi-squared_distribution</a>.
922      *
923      * @param {Number} k <i>> 0</i>
924      * @returns Number
925      * @memberof JXG.Math.Statistics
926      */
927     randomChisquare: function (nu) {
928         // Knuth, TAOCP 2, p 130
929 
930         if (nu <= 0) {
931             return NaN;
932         }
933 
934         return 2 * this.randomGamma(nu * 0.5);
935     },
936 
937     /**
938      * Generate value of a random variable with F-distribution with d<sub>1</sub> and d<sub>2</sub> degrees of freedom.
939      * See <a href="https://en.wikipedia.org/wiki/F-distribution">https://en.wikipedia.org/wiki/F-distribution</a>.
940      * @param {Number} d1 <i>> 0</i>
941      * @param {Number} d2 <i>> 0</i>
942      * @returns Number
943      * @memberof JXG.Math.Statistics
944      */
945     randomF: function (nu1, nu2) {
946         // Knuth, TAOCP 2, p 130
947         var y1, y2;
948 
949         if (nu1 <= 0 || nu2 <= 0) {
950             return NaN;
951         }
952 
953         y1 = this.randomChisquare(nu1);
954         y2 = this.randomChisquare(nu2);
955 
956         return (y1 * nu2) / (y2 * nu1);
957     },
958 
959     /**
960      * Generate value of a random variable with Students-t-distribution with ν degrees of freedom.
961      * See <a href="https://en.wikipedia.org/wiki/Student%27s_t-distribution">https://en.wikipedia.org/wiki/Student%27s_t-distribution</a>.
962      * @param {Number} nu <i>> 0</i>
963      * @returns Number
964      * @memberof JXG.Math.Statistics
965      */
966     randomT: function (nu) {
967         // Knuth, TAOCP 2, p 130
968         var y1, y2;
969 
970         if (nu <= 0) {
971             return NaN;
972         }
973 
974         y1 = this.randomNormal(0, 1);
975         y2 = this.randomChisquare(nu);
976 
977         return y1 / Math.sqrt(y2 / nu);
978     },
979 
980     /**
981      * Generate values for a random variable in binomial distribution with parameters <i>n</i> and <i>p</i>.
982      * See <a href="https://en.wikipedia.org/wiki/Binomial_distribution">https://en.wikipedia.org/wiki/Binomial_distribution</a>.
983      * It uses algorithm BG from <a href="https://dl.acm.org/doi/pdf/10.1145/42372.42381">https://dl.acm.org/doi/pdf/10.1145/42372.42381</a>.
984      *
985      * @param {Number} n Number of trials (n >= 0)
986      * @param {Number} p Probability (0 <= p <= 1)
987      * @returns Number Integer value of a random variable in binomial distribution
988      * @memberof JXG.Math.Statistics
989      *
990      * @example
991      * let board = JXG.JSXGraph.initBoard('jxgbox',
992      *     { boundingbox: [-1.7, .5, 30, -.03], axis: true });
993      *
994      * let runs = [
995      *     [0.5, 20, 'blue'],
996      *     [0.7, 20, 'green'],
997      *     [0.5, 40, 'red'],
998      * ];
999      *
1000      * let labelY = .4;
1001      * runs.forEach((run, i) => {
1002      *     board.create('segment', [[7, labelY - (i / 50)], [9, labelY - (i / 50)]], { strokeColor: run[2] });
1003      *     board.create('text', [10, labelY - (i / 50), `p=${run[0]}, n=${run[1]}`]);
1004      *
1005      *     let x = Array(50000).fill(0).map(() => JXG.Math.Statistics.randomBinomial(run[1], run[0]));
1006      *     let res = JXG.Math.Statistics.histogram(x, {
1007      *         bins: 40,
1008      *         density: true,
1009      *         cumulative: false,
1010      *         range: [0, 40]
1011      *     });
1012      *     board.create('curve', [res[1], res[0]], { strokeColor: run[2] });
1013      * });
1014      *
1015      *
1016      * </pre><div id="JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302-3" class="jxgbox" style="width: 300px; height: 300px;"></div>
1017      * <script type="text/javascript">
1018      * {
1019      *  let board = JXG.JSXGraph.initBoard('JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302-3',
1020      *       { boundingbox: [-1.7, .5, 30, -.03], axis: true});
1021      *
1022      *   let runs = [
1023      *       [0.5, 20, 'blue'],
1024      *       [0.7, 20, 'green'],
1025      *       [0.5, 40, 'red'],
1026      *   ]
1027      *
1028      *   let labelY = .4
1029      *   runs.forEach((run,i) => {
1030      *       board.create('segment',[[7,labelY-(i/50)],[9,labelY-(i/50)]],{strokeColor:run[2]})
1031      *       board.create('text',[10,labelY-(i/50),`p=${run[0]}, n=${run[1]}`])
1032      *
1033      *       let x = Array(50000).fill(0).map(() => JXG.Math.Statistics.randomBinomial(run[1],run[0]))
1034      *       let res = JXG.Math.Statistics.histogram(x, { bins: 40, density: true, cumulative: false, range: [0, 40] });
1035      *       board.create('curve', [res[1], res[0]], { strokeColor: run[2]});
1036      *   })
1037      * }
1038      * </script><pre>
1039      *
1040      */
1041     randomBinomial: function (n, p) {
1042         var x, y, c,
1043             a, b, N1;
1044 
1045         if (p < 0 || p > 1 || n < 0) {
1046             return NaN;
1047         }
1048 
1049         // Edge cases
1050         if (p === 0) {
1051             return 0;
1052         }
1053         if (p === 1) {
1054             return n;
1055         }
1056 
1057         // Now, we can assume 0 < p < 1.
1058 
1059         // Fast path for common cases
1060         if (n === 0) {
1061             return 0;
1062         }
1063         if (n === 1) {
1064             return ((Math.random() < p) ? 1 : 0);
1065         }
1066 
1067         // Exploit symmetry
1068         if (p > 0.5) {
1069             return n - this.randomBinomial(n, 1 - p);
1070         }
1071 
1072         // General case: n > 1, p <= 0.5
1073         if (n < 100) {
1074             // n small:
1075             // Algorithm BG (Devroye) from:
1076             // https://dl.acm.org/doi/pdf/10.1145/42372.42381
1077             // Time O(np) so suitable for np small only.
1078             x = -1;
1079             y = 0;
1080 
1081             c = Math.log(1 - p);
1082             if (c === 0) {
1083                 return 0;
1084             }
1085 
1086             do {
1087                 x += 1;
1088                 y += Math.floor(Math.log(Math.random()) / c) + 1;
1089             } while (y < n);
1090         } else {
1091             // n large:
1092             // Knuth, TAOCP 2, p 131
1093             a = 1 + Math.floor(n * 0.5);
1094             b = n - a + 1;
1095             x = this.randomBeta(a, b);
1096             if (x >= p) {
1097                 N1 = this.randomBinomial(a - 1, p / x);
1098                 x = N1;
1099             } else {
1100                 N1 = this.randomBinomial(b - 1, (p - x) / (1 - x));
1101                 x = a + N1;
1102             }
1103         }
1104         return x;
1105     },
1106 
1107     /**
1108      * Generate values for a random variable in geometric distribution with probability <i>p</i>.
1109      * See <a href="https://en.wikipedia.org/wiki/Geometric_distribution">https://en.wikipedia.org/wiki/Geometric_distribution</a>.
1110      *
1111      * @param {Number} p (0 <= p <= 1)
1112      * @returns Number
1113      * @memberof JXG.Math.Statistics
1114      */
1115     randomGeometric: function (p) {
1116         var u;
1117 
1118         if (p < 0 || p > 1) {
1119             return NaN;
1120         }
1121         // Knuth, TAOCP 2, p 131
1122         u = Math.random();
1123 
1124         return Math.ceil(Math.log(u) / Math.log(1 - p));
1125     },
1126 
1127     /**
1128      * Generate values for a random variable in Poisson distribution with mean <i>mu</i>.
1129      * See <a href="https://en.wikipedia.org/wiki/Poisson_distribution">https://en.wikipedia.org/wiki/Poisson_distribution</a>.
1130      *
1131      * @param {Number} mu (0 < mu)
1132      * @returns Number
1133      * @memberof JXG.Math.Statistics
1134      */
1135     randomPoisson: function (mu) {
1136         var e = Math.exp(-mu),
1137             N,
1138             m = 0,
1139             u = 1,
1140             x,
1141             alpha = 7 / 8;
1142 
1143         if (mu <= 0) {
1144             return NaN;
1145         }
1146 
1147         // Knuth, TAOCP 2, p 132
1148         if (mu < 10) {
1149             do {
1150                 u *= Math.random();
1151                 m += 1;
1152             } while (u > e);
1153             N = m - 1;
1154         } else {
1155             m = Math.floor(alpha * mu);
1156             x = this.randomGamma(m);
1157             if (x < mu) {
1158                 N = m + this.randomPoisson(mu - x);
1159             } else {
1160                 N = this.randomBinomial(m - 1, mu / x);
1161             }
1162         }
1163         return N;
1164     },
1165 
1166     /**
1167      * Generate values for a random variable in Pareto distribution with
1168      * shape <i>gamma</i> and scale <i>k</i>.
1169      * See <a href="https://en.wikipedia.org/wiki/Pareto_distribution">https://en.wikipedia.org/wiki/Pareto_distribution</a>.
1170      * Method: use inverse transformation sampling.
1171      *
1172      * @param {Number} gamma shape (0 < gamma)
1173      * @param {Number} k scale (0 < k < x)
1174      * @returns Number
1175      * @memberof JXG.Math.Statistics
1176      */
1177     randomPareto: function (gamma, k) {
1178         var u = Math.random();
1179 
1180         if (gamma <= 0 || k <= 0) {
1181             return NaN;
1182         }
1183         return k * Math.pow(1 - u, -1 / gamma);
1184     },
1185 
1186     /**
1187      * Generate values for a random variable in hypergeometric distribution.
1188      * Samples are drawn from a hypergeometric distribution with specified parameters, <i>good</i> (ways to make a good selection),
1189      * <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>).
1190      * <p>
1191      * Naive implementation with runtime <i>O(samples)</i>.
1192      *
1193      * @param {Number} good ways to make a good selection
1194      * @param {Number} bad ways to make a bad selection
1195      * @param {Number} samples number of items sampled
1196      * @returns
1197      * @memberof JXG.Math.Statistics
1198      */
1199     randomHypergeometric: function (good, bad, k) {
1200         var i, u,
1201             x = 0,
1202             // kk,
1203             // n = good + bad,
1204             d1 = good + bad - k,
1205             d2 = Math.min(good, bad),
1206             y = d2;
1207 
1208         if (good < 1 || bad < 1 || k > good + bad) {
1209             return NaN;
1210         }
1211 
1212         // Naive method
1213         // kk = Math.min(k, n - k);
1214         // for (i = 0; i < k; i ++) {
1215         //     u = Math.random();
1216         //     if (n * u <= good) {
1217         //         x += 1;
1218         //         if (x === good) {
1219         //             return x;
1220         //         }
1221         //         good -= 1;
1222         //     }
1223         //     n -= 1;
1224         // }
1225         // return x;
1226 
1227         // Implementation from
1228         // Monte Carlo by George S. Fishman
1229         // https://link.springer.com/book/10.1007/978-1-4757-2553-7
1230         // page 218
1231         //
1232         i = k;
1233         while (y * i > 0) {
1234             u = Math.random();
1235             y -= Math.floor(u + y / (d1 + i));
1236             i -= 1;
1237         }
1238         x = d2 - y;
1239         if (good <= bad) {
1240             return x;
1241         } else {
1242             return k - x;
1243         }
1244     },
1245 
1246     /**
1247      * Compute the histogram of a dataset.
1248      * Optional parameters can be supplied through a JavaScript object
1249      * with the following default values:
1250      * <pre>
1251      * {
1252      *   bins: 10,          // Number of bins
1253      *   range: false,      // false or array. The lower and upper range of the bins.
1254      *                      // If not provided, range is simply [min(x), max(x)].
1255      *                      // Values outside the range are ignored.
1256      *   density: false,    // If true, normalize the counts by dividing by sum(counts)
1257      *   cumulative: false
1258      * }
1259      * </pre>
1260      * The function returns an array containing two arrays. The first array is of length bins+1
1261      * containing the start values of the bins. The last entry contains the end values of the last bin.
1262      * <p>
1263      * The second array contains the counts of each bin.
1264      * @param {Array} x
1265      * @param {Object} opt Optional parameters
1266      * @returns Array [bin, counts] Array bins contains start values of bins, array counts contains
1267      * the number of entries of x which are contained in each bin.
1268      * @memberof JXG.Math.Statistics
1269      *
1270      * @example
1271      *  let board = JXG.JSXGraph.initBoard('jxgbox',
1272      *       { boundingbox: [-1.7, .5, 20, -.03], axis: true});
1273      *  let board2 = JXG.JSXGraph.initBoard('jxgbox2',
1274      *       { boundingbox: [-1.6, 1.1, 20, -.06], axis: true});
1275      *
1276      *   let runs = [
1277      *       [0.5, 1.0, 'brown'],
1278      *       [1.0, 2.0, 'red'],
1279      *       [2.0, 2.0, 'orange'],
1280      *       [3.0, 2.0, 'yellow'],
1281      *       [5.0, 1.0, 'green'],
1282      *       [9.0, 0.5, 'black'],
1283      *       [7.5, 1.0, 'purple'],
1284      *   ]
1285      *
1286      *   let labelY = .4
1287      *   runs.forEach((run,i) => {
1288      *       board.create('segment',[[7,labelY-(i/50)],[9,labelY-(i/50)]],{strokeColor:run[2]})
1289      *       board.create('text',[10,labelY-(i/50),`k=${run[0]}, θ=${run[1]}`])
1290      *
1291      *       // density
1292      *       let x = Array(50000).fill(0).map(() => JXG.Math.Statistics.randomGamma(run[0],run[1]))
1293      *       let res = JXG.Math.Statistics.histogram(x, { bins: 50, density: true, cumulative: false, range: [0, 20] });
1294      *       board.create('curve', [res[1], res[0]], { strokeColor: run[2], strokeWidth:2});
1295      *
1296      *       // cumulative density
1297      *       res = JXG.Math.Statistics.histogram(x, { bins: 50, density: true, cumulative: true, range: [0, 20] });
1298      *       res[0].unshift(0)  // add zero to front so cumulative starts at zero
1299      *       res[1].unshift(0)
1300      *       board2.create('curve', [res[1], res[0]], { strokeColor: run[2], strokeWidth:2 });
1301      *   })
1302      *
1303      *
1304      * </pre><div id="JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302" class="jxgbox" style="width: 300px; height: 300px; float:left;"></div>
1305      * <div style='float:left;'>  </div>
1306      * <div id="JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302-2" class="jxgbox" style="width: 300px; height: 300px;"></div>
1307      * <script type="text/javascript">
1308      * {
1309      *  let board = JXG.JSXGraph.initBoard('JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302',
1310      *       { boundingbox: [-1.7, .5, 20, -.03], axis: true});
1311      *  let board2 = JXG.JSXGraph.initBoard('JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302-2',
1312      *       { boundingbox: [-1.6, 1.1, 20, -.06], axis: true});
1313      *
1314      *   let runs = [
1315      *       [0.5, 1.0, 'brown'],
1316      *       [1.0, 2.0, 'red'],
1317      *       [2.0, 2.0, 'orange'],
1318      *       [3.0, 2.0, 'yellow'],
1319      *       [5.0, 1.0, 'green'],
1320      *       [9.0, 0.5, 'black'],
1321      *       [7.5, 1.0, 'purple'],
1322      *   ]
1323      *
1324      *   let labelY = .4
1325      *   runs.forEach((run,i) => {
1326      *       board.create('segment',[[7,labelY-(i/50)],[9,labelY-(i/50)]],{strokeColor:run[2]})
1327      *       board.create('text',[10,labelY-(i/50),`k=${run[0]}, θ=${run[1]}`])
1328      *
1329      *       let x = Array(50000).fill(0).map(() => JXG.Math.Statistics.randomGamma(run[0],run[1]))
1330      *       let res = JXG.Math.Statistics.histogram(x, { bins: 50, density: true, cumulative: false, range: [0, 20] });
1331      *       board.create('curve', [res[1], res[0]], { strokeColor: run[2], strokeWidth:2});
1332      *
1333      *       // cumulative density
1334      *       res = JXG.Math.Statistics.histogram(x, { bins: 50, density: true, cumulative: true, range: [0, 20] });
1335      *       res[0].unshift(0)  // add zero to front so cumulative starts at zero
1336      *       res[1].unshift(0)
1337      *       board2.create('curve', [res[1], res[0]], { strokeColor: run[2], strokeWidth:2 });
1338      *   })
1339      * }
1340      * </script><pre>
1341      *
1342      */
1343     histogram: function (x, opt) {
1344         var i, le, k,
1345             mi, ma, num_bins, delta,
1346             range,
1347             s,
1348             counts = [],
1349             bins = [],
1350             no_bin = 0;        // Count of long tail elements not in histogram range
1351 
1352         // Evaluate number of bins
1353         num_bins = opt.bins || 10;
1354 
1355         // Evaluate range
1356         range = opt.range || false;
1357         if (range === false) {
1358             mi = Math.min.apply(null, x);
1359             ma = Math.max.apply(null, x);
1360         } else {
1361             mi = range[0];
1362             ma = range[1];
1363         }
1364 
1365         // Set uniform delta
1366         if (num_bins > 0) {
1367             delta = (ma - mi) / (num_bins - 1);
1368         } else {
1369             delta = 0;
1370         }
1371 
1372         // Set the bins and init the counts array
1373         for (i = 0; i < num_bins; i++) {
1374             counts.push(0);
1375             bins.push(mi + i * delta);
1376         }
1377         // bins.push(ma);
1378 
1379         // Determine the counts
1380         le = x.length;
1381         for (i = 0; i < le; i++) {
1382             k = Math.floor((x[i] - mi) / delta);
1383             if (k >= 0 && k < num_bins) {
1384                 counts[k] += 1;
1385             } else {
1386                 no_bin += 1;
1387             }
1388         }
1389 
1390         // Normalize if density===true
1391         if (opt.density) {
1392             s = JXG.Math.Statistics.sum(counts) + no_bin; // Normalize including long tail
1393             for (i = 0; i < num_bins; i++) {
1394                 counts[i] /= (s * delta);
1395                 // counts[i] /= s;
1396             }
1397         }
1398 
1399         // Cumulative counts
1400         if (opt.cumulative) {
1401             if (opt.density) {
1402                 for (i = 0; i < num_bins; i++) {
1403                     counts[i] *= delta;  // Normalize
1404                 }
1405             } for (i = 1; i < num_bins; i++) {
1406                 counts[i] += counts[i - 1];
1407             }
1408         }
1409 
1410         return [counts, bins];
1411     }
1412 };
1413 
1414 export default Mat.Statistics;
1415