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, { bins: 40, density: true, cumulative: false, range: [0, 40] });
1007      *       board.create('curve', [res[1], res[0]], { strokeColor: run[2]});
1008      *   })
1009      * }
1010      *
1011      *
1012      * </pre><div id="JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302-3" class="jxgbox" style="width: 300px; height: 300px;"></div>
1013      * <script type="text/javascript">
1014      * {
1015      *  let board = JXG.JSXGraph.initBoard('JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302-3',
1016      *       { boundingbox: [-1.7, .5, 30, -.03], axis: true});
1017      *
1018      *   let runs = [
1019      *       [0.5, 20, 'blue'],
1020      *       [0.7, 20, 'green'],
1021      *       [0.5, 40, 'red'],
1022      *   ]
1023      *
1024      *   let labelY = .4
1025      *   runs.forEach((run,i) => {
1026      *       board.create('segment',[[7,labelY-(i/50)],[9,labelY-(i/50)]],{strokeColor:run[2]})
1027      *       board.create('text',[10,labelY-(i/50),`p=${run[0]}, n=${run[1]}`])
1028      *
1029      *       let x = Array(50000).fill(0).map(() => JXG.Math.Statistics.randomBinomial(run[1],run[0]))
1030      *       let res = JXG.Math.Statistics.histogram(x, { bins: 40, density: true, cumulative: false, range: [0, 40] });
1031      *       board.create('curve', [res[1], res[0]], { strokeColor: run[2]});
1032      *   })
1033      * }
1034      * </script><pre>
1035      *
1036      */
1037     randomBinomial: function (n, p) {
1038         var x, y, c,
1039             a, b, N1;
1040 
1041         if (p < 0 || p > 1 || n < 0) {
1042             return NaN;
1043         }
1044 
1045         // Edge cases
1046         if (p === 0) {
1047             return 0;
1048         }
1049         if (p === 1) {
1050             return n;
1051         }
1052 
1053         // Now, we can assume 0 < p < 1.
1054 
1055         // Fast path for common cases
1056         if (n === 0) {
1057             return 0;
1058         }
1059         if (n === 1) {
1060             return ((Math.random() < p) ? 1 : 0);
1061         }
1062 
1063         // Exploit symmetry
1064         if (p > 0.5) {
1065             return n - this.randomBinomial(n, 1 - p);
1066         }
1067 
1068         // General case: n > 1, p <= 0.5
1069         if (n < 100) {
1070             // n small:
1071             // Algorithm BG (Devroye) from:
1072             // https://dl.acm.org/doi/pdf/10.1145/42372.42381
1073             // Time O(np) so suitable for np small only.
1074             x = -1;
1075             y = 0;
1076 
1077             c = Math.log(1 - p);
1078             if (c === 0) {
1079                 return 0;
1080             }
1081 
1082             do {
1083                 x += 1;
1084                 y += Math.floor(Math.log(Math.random()) / c) + 1;
1085             } while (y < n);
1086         } else {
1087             // n large:
1088             // Knuth, TAOCP 2, p 131
1089             a = 1 + Math.floor(n * 0.5);
1090             b = n - a + 1;
1091             x = this.randomBeta(a, b);
1092             if (x >= p) {
1093                 N1 = this.randomBinomial(a - 1, p / x);
1094                 x = N1;
1095             } else {
1096                 N1 = this.randomBinomial(b - 1, (p - x) / (1 - x));
1097                 x = a + N1;
1098             }
1099         }
1100         return x;
1101     },
1102 
1103     /**
1104      * Generate values for a random variable in geometric distribution with probability <i>p</i>.
1105      * See <a href="https://en.wikipedia.org/wiki/Geometric_distribution">https://en.wikipedia.org/wiki/Geometric_distribution</a>.
1106      *
1107      * @param {Number} p (0 <= p <= 1)
1108      * @returns Number
1109      * @memberof JXG.Math.Statistics
1110      */
1111     randomGeometric: function (p) {
1112         var u;
1113 
1114         if (p < 0 || p > 1) {
1115             return NaN;
1116         }
1117         // Knuth, TAOCP 2, p 131
1118         u = Math.random();
1119 
1120         return Math.ceil(Math.log(u) / Math.log(1 - p));
1121     },
1122 
1123     /**
1124      * Generate values for a random variable in Poisson distribution with mean <i>mu</i>.
1125      * See <a href="https://en.wikipedia.org/wiki/Poisson_distribution">https://en.wikipedia.org/wiki/Poisson_distribution</a>.
1126      *
1127      * @param {Number} mu (0 < mu)
1128      * @returns Number
1129      * @memberof JXG.Math.Statistics
1130      */
1131     randomPoisson: function (mu) {
1132         var e = Math.exp(-mu),
1133             N,
1134             m = 0,
1135             u = 1,
1136             x,
1137             alpha = 7 / 8;
1138 
1139         if (mu <= 0) {
1140             return NaN;
1141         }
1142 
1143         // Knuth, TAOCP 2, p 132
1144         if (mu < 10) {
1145             do {
1146                 u *= Math.random();
1147                 m += 1;
1148             } while (u > e);
1149             N = m - 1;
1150         } else {
1151             m = Math.floor(alpha * mu);
1152             x = this.randomGamma(m);
1153             if (x < mu) {
1154                 N = m + this.randomPoisson(mu - x);
1155             } else {
1156                 N = this.randomBinomial(m - 1, mu / x);
1157             }
1158         }
1159         return N;
1160     },
1161 
1162     /**
1163      * Generate values for a random variable in Pareto distribution with
1164      * shape <i>gamma</i> and scale <i>k</i>.
1165      * See <a href="https://en.wikipedia.org/wiki/Pareto_distribution">https://en.wikipedia.org/wiki/Pareto_distribution</a>.
1166      * Method: use inverse transformation sampling.
1167      *
1168      * @param {Number} gamma shape (0 < gamma)
1169      * @param {Number} k scale (0 < k < x)
1170      * @returns Number
1171      * @memberof JXG.Math.Statistics
1172      */
1173     randomPareto: function (gamma, k) {
1174         var u = Math.random();
1175 
1176         if (gamma <= 0 || k <= 0) {
1177             return NaN;
1178         }
1179         return k * Math.pow(1 - u, -1 / gamma);
1180     },
1181 
1182     /**
1183      * Generate values for a random variable in hypergeometric distribution.
1184      * Samples are drawn from a hypergeometric distribution with specified parameters, <i>good</i> (ways to make a good selection),
1185      * <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>).
1186      * <p>
1187      * Naive implementation with runtime <i>O(samples)</i>.
1188      *
1189      * @param {Number} good ways to make a good selection
1190      * @param {Number} bad ways to make a bad selection
1191      * @param {Number} samples number of items sampled
1192      * @returns
1193      * @memberof JXG.Math.Statistics
1194      */
1195     randomHypergeometric: function (good, bad, k) {
1196         var i, u,
1197             x = 0,
1198             // kk,
1199             // n = good + bad,
1200             d1 = good + bad - k,
1201             d2 = Math.min(good, bad),
1202             y = d2;
1203 
1204         if (good < 1 || bad < 1 || k > good + bad) {
1205             return NaN;
1206         }
1207 
1208         // Naive method
1209         // kk = Math.min(k, n - k);
1210         // for (i = 0; i < k; i ++) {
1211         //     u = Math.random();
1212         //     if (n * u <= good) {
1213         //         x += 1;
1214         //         if (x === good) {
1215         //             return x;
1216         //         }
1217         //         good -= 1;
1218         //     }
1219         //     n -= 1;
1220         // }
1221         // return x;
1222 
1223         // Implementation from
1224         // Monte Carlo by George S. Fishman
1225         // https://link.springer.com/book/10.1007/978-1-4757-2553-7
1226         // page 218
1227         //
1228         i = k;
1229         while (y * i > 0) {
1230             u = Math.random();
1231             y -= Math.floor(u + y / (d1 + i));
1232             i -= 1;
1233         }
1234         x = d2 - y;
1235         if (good <= bad) {
1236             return x;
1237         } else {
1238             return k - x;
1239         }
1240     },
1241 
1242     /**
1243      * Compute the histogram of a dataset.
1244      * Optional parameters can be supplied through a JavaScript object
1245      * with the following default values:
1246      * <pre>
1247      * {
1248      *   bins: 10,          // Number of bins
1249      *   range: false,      // false or array. The lower and upper range of the bins.
1250      *                      // If not provided, range is simply [min(x), max(x)].
1251      *                      // Values outside the range are ignored.
1252      *   density: false,    // If true, normalize the counts by dividing by sum(counts)
1253      *   cumulative: false
1254      * }
1255      * </pre>
1256      * The function returns an array containing two arrays. The first array is of length bins+1
1257      * containing the start values of the bins. The last entry contains the end values of the last bin.
1258      * <p>
1259      * The second array contains the counts of each bin.
1260      * @param {Array} x
1261      * @param {Object} opt Optional parameters
1262      * @returns Array [bin, counts] Array bins contains start values of bins, array counts contains
1263      * the number of entries of x which are contained in each bin.
1264      * @memberof JXG.Math.Statistics
1265      *
1266      * @example
1267      *  let board = JXG.JSXGraph.initBoard('jxgbox',
1268      *       { boundingbox: [-1.7, .5, 20, -.03], axis: true});
1269      *  let board2 = JXG.JSXGraph.initBoard('jxgbox2',
1270      *       { boundingbox: [-1.6, 1.1, 20, -.06], axis: true});
1271      *
1272      *   let runs = [
1273      *       [0.5, 1.0, 'brown'],
1274      *       [1.0, 2.0, 'red'],
1275      *       [2.0, 2.0, 'orange'],
1276      *       [3.0, 2.0, 'yellow'],
1277      *       [5.0, 1.0, 'green'],
1278      *       [9.0, 0.5, 'black'],
1279      *       [7.5, 1.0, 'purple'],
1280      *   ]
1281      *
1282      *   let labelY = .4
1283      *   runs.forEach((run,i) => {
1284      *       board.create('segment',[[7,labelY-(i/50)],[9,labelY-(i/50)]],{strokeColor:run[2]})
1285      *       board.create('text',[10,labelY-(i/50),`k=${run[0]}, θ=${run[1]}`])
1286      *
1287      *       // density
1288      *       let x = Array(50000).fill(0).map(() => JXG.Math.Statistics.randomGamma(run[0],run[1]))
1289      *       let res = JXG.Math.Statistics.histogram(x, { bins: 50, density: true, cumulative: false, range: [0, 20] });
1290      *       board.create('curve', [res[1], res[0]], { strokeColor: run[2], strokeWidth:2});
1291      *
1292      *       // cumulative density
1293      *       res = JXG.Math.Statistics.histogram(x, { bins: 50, density: true, cumulative: true, range: [0, 20] });
1294      *       res[0].unshift(0)  // add zero to front so cumulative starts at zero
1295      *       res[1].unshift(0)
1296      *       board2.create('curve', [res[1], res[0]], { strokeColor: run[2], strokeWidth:2 });
1297      *   })
1298      *
1299      *
1300      * </pre><div id="JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302" class="jxgbox" style="width: 300px; height: 300px; float:left;"></div>
1301      * <div style='float:left;'>  </div>
1302      * <div id="JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302-2" class="jxgbox" style="width: 300px; height: 300px;"></div>
1303      * <script type="text/javascript">
1304      * {
1305      *  let board = JXG.JSXGraph.initBoard('JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302',
1306      *       { boundingbox: [-1.7, .5, 20, -.03], axis: true});
1307      *  let board2 = JXG.JSXGraph.initBoard('JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302-2',
1308      *       { boundingbox: [-1.6, 1.1, 20, -.06], axis: true});
1309      *
1310      *   let runs = [
1311      *       [0.5, 1.0, 'brown'],
1312      *       [1.0, 2.0, 'red'],
1313      *       [2.0, 2.0, 'orange'],
1314      *       [3.0, 2.0, 'yellow'],
1315      *       [5.0, 1.0, 'green'],
1316      *       [9.0, 0.5, 'black'],
1317      *       [7.5, 1.0, 'purple'],
1318      *   ]
1319      *
1320      *   let labelY = .4
1321      *   runs.forEach((run,i) => {
1322      *       board.create('segment',[[7,labelY-(i/50)],[9,labelY-(i/50)]],{strokeColor:run[2]})
1323      *       board.create('text',[10,labelY-(i/50),`k=${run[0]}, θ=${run[1]}`])
1324      *
1325      *       let x = Array(50000).fill(0).map(() => JXG.Math.Statistics.randomGamma(run[0],run[1]))
1326      *       let res = JXG.Math.Statistics.histogram(x, { bins: 50, density: true, cumulative: false, range: [0, 20] });
1327      *       board.create('curve', [res[1], res[0]], { strokeColor: run[2], strokeWidth:2});
1328      *
1329      *       // cumulative density
1330      *       res = JXG.Math.Statistics.histogram(x, { bins: 50, density: true, cumulative: true, range: [0, 20] });
1331      *       res[0].unshift(0)  // add zero to front so cumulative starts at zero
1332      *       res[1].unshift(0)
1333      *       board2.create('curve', [res[1], res[0]], { strokeColor: run[2], strokeWidth:2 });
1334      *   })
1335      * }
1336      * </script><pre>
1337      *
1338      */
1339     histogram: function (x, opt) {
1340         var i, le, k,
1341             mi, ma, num_bins, delta,
1342             range,
1343             s,
1344             counts = [],
1345             bins = [];
1346 
1347         // Evaluate number of bins
1348         num_bins = opt.bins || 10;
1349 
1350         // Evaluate range
1351         range = opt.range || false;
1352         if (range === false) {
1353             mi = Math.min.apply(null, x);
1354             ma = Math.max.apply(null, x);
1355         } else {
1356             mi = range[0];
1357             ma = range[1];
1358         }
1359 
1360         // Set uniform delta
1361         if (num_bins > 0) {
1362             delta = (ma - mi) / (num_bins - 1);
1363         } else {
1364             delta = 0;
1365         }
1366 
1367         // Set the bins and init the counts array
1368         for (i = 0; i < num_bins; i++) {
1369             counts.push(0);
1370             bins.push(mi + i * delta);
1371         }
1372         // bins.push(ma);
1373 
1374         // Determine the counts
1375         le = x.length;
1376         for (i = 0; i < le; i++) {
1377             k = Math.floor((x[i] - mi) / delta);
1378             if (k >= 0 && k < num_bins) {
1379                 counts[k] += 1;
1380             }
1381         }
1382 
1383         // Normalize if density===true
1384         if (opt.density) {
1385             s = JXG.Math.Statistics.sum(counts);
1386             for (i = 0; i < num_bins; i++) {
1387                 counts[i] /= (s * delta);
1388                 // counts[i] /= s;
1389             }
1390         }
1391 
1392         // Cumulative counts
1393         if (opt.cumulative) {
1394             if (opt.density) {
1395                 for (i = 0; i < num_bins; i++) {
1396                     counts[i] *= delta;  // normalize
1397                 }
1398             } for (i = 1; i < num_bins; i++) {
1399                 counts[i] += counts[i - 1];
1400             }
1401         }
1402 
1403         return [counts, bins];
1404     }
1405 };
1406 
1407 export default Mat.Statistics;
1408