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