1 /*
  2     Copyright 2008-2022
  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 <http://www.gnu.org/licenses/>
 29     and <http://opensource.org/licenses/MIT/>.
 30  */
 31 
 32 
 33 /*global JXG: true, define: true*/
 34 /*jslint nomen: true, plusplus: true*/
 35 
 36 /* depends:
 37  jxg
 38  math/math
 39  utils/type
 40  */
 41 
 42 define(['jxg', 'math/math', 'utils/type'], function (JXG, Mat, Type) {
 43 
 44     "use strict";
 45 
 46     /**
 47      * Functions for mathematical statistics. Most functions are like in the statistics package R.
 48      * @name JXG.Math.Statistics
 49      * @exports Mat.Statistics as JXG.Math.Statistics
 50      * @namespace
 51      */
 52     Mat.Statistics = {
 53         /**
 54          * Sums up all elements of the given array.
 55          * @param {Array} arr An array of numbers.
 56          * @returns {Number}
 57          * @memberof JXG.Math.Statistics
 58          */
 59         sum: function (arr) {
 60             var i,
 61                 len = arr.length,
 62                 res = 0;
 63 
 64             for (i = 0; i < len; i++) {
 65                 res += arr[i];
 66             }
 67             return res;
 68         },
 69 
 70         /**
 71          * Multiplies all elements of the given array.
 72          * @param {Array} arr An array of numbers.
 73          * @returns {Number}
 74          * @memberof JXG.Math.Statistics
 75          */
 76         prod: function (arr) {
 77             var i,
 78                 len = arr.length,
 79                 res = 1;
 80 
 81             for (i = 0; i < len; i++) {
 82                 res *= arr[i];
 83             }
 84             return res;
 85         },
 86 
 87         /**
 88          * Determines the mean value of the values given in an array.
 89          * @param {Array} arr
 90          * @returns {Number}
 91          * @memberof JXG.Math.Statistics
 92          */
 93         mean: function (arr) {
 94             if (arr.length > 0) {
 95                 return this.sum(arr) / arr.length;
 96             }
 97 
 98             return 0.0;
 99         },
100 
101         /**
102          * The median of a finite set of values is the value that divides the set
103          * into two equal sized subsets.
104          * @param {Array} arr The set of values.
105          * @returns {Number}
106          * @memberof JXG.Math.Statistics
107          */
108         median: function (arr) {
109             var tmp, len;
110 
111             if (arr.length > 0) {
112                 if (ArrayBuffer.isView(arr)) {
113                     tmp = new Float64Array(arr);
114                     tmp.sort();
115                 } else {
116                     tmp = arr.slice(0);
117                     tmp.sort(function (a, b) {
118                         return a - b;
119                     });
120                 }
121                 len = tmp.length;
122 
123                 if (len & 1) { // odd
124                     return tmp[parseInt(len * 0.5, 10)];
125                 }
126 
127                 return (tmp[len * 0.5 - 1] + tmp[len * 0.5]) * 0.5;
128             }
129 
130             return 0.0;
131         },
132 
133         /**
134          * The P-th percentile ( 0 < P ≤ 100 ) of a list of N ordered values (sorted from least to greatest)
135          * is the smallest value in the list such that no more than P percent of the data is strictly less
136          * than the value and at least P percent of the data is less than or equal to that value. See {@link https://en.wikipedia.org/wiki/Percentile}.
137          *
138          * Here, the <i>linear interpolation between closest ranks</i> method is used.
139          * @param {Array} arr The set of values, need not be ordered.
140          * @param {Number|Array} percentile One or several percentiles
141          * @returns {Number|Array} Depending if a number or an array is the input for percentile, a number or an array containing the percentils
142          * is returned.
143          */
144         percentile: function(arr, percentile) {
145             var tmp, len, i, p, res = [], 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, res, i, len = arr.length;
193 
194             if (len > 1) {
195                 m = this.mean(arr);
196                 res = 0;
197                 for (i = 0; i < len; i++) {
198                     res += (arr[i] - m) * (arr[i] - m);
199                 }
200                 return res / (arr.length - 1);
201             }
202 
203             return 0.0;
204         },
205 
206         /**
207          * Determines the <strong>s</strong>tandard <strong>d</strong>eviation which shows how much
208          * variation there is from the average value of a set of numbers.
209          * @param {Array} arr
210          * @returns {Number}
211          * @memberof JXG.Math.Statistics
212          */
213         sd: function (arr) {
214             return Math.sqrt(this.variance(arr));
215         },
216 
217         /**
218          * Weighted mean value is basically the same as {@link JXG.Math.Statistics.mean} but here the values
219          * are weighted, i.e. multiplied with another value called <em>weight</em>. The weight values are given
220          * as a second array with the same length as the value array..
221          * @throws {Error} If the dimensions of the arrays don't match.
222          * @param {Array} arr Set of alues.
223          * @param {Array} w Weight values.
224          * @returns {Number}
225          * @memberof JXG.Math.Statistics
226          */
227         weightedMean: function (arr, w) {
228             if (arr.length !== w.length) {
229                 throw new Error('JSXGraph error (Math.Statistics.weightedMean): Array dimension mismatch.');
230             }
231 
232             if (arr.length > 0) {
233                 return this.mean(this.multiply(arr, w));
234             }
235 
236             return 0.0;
237         },
238 
239         /**
240          * Extracts the maximum value from the array.
241          * @param {Array} arr
242          * @returns {Number} The highest number from the array. It returns <tt>NaN</tt> if not every element could be
243          * interpreted as a number and <tt>-Infinity</tt> if an empty array is given or no element could be interpreted
244          * as a number.
245          * @memberof JXG.Math.Statistics
246          */
247         max: function (arr) {
248             return Math.max.apply(this, arr);
249         },
250 
251         /**
252          * Extracts the minimum value from the array.
253          * @param {Array} arr
254          * @returns {Number} The lowest number from the array. It returns <tt>NaN</tt> if not every element could be
255          * interpreted as a number and <tt>Infinity</tt> if an empty array is given or no element could be interpreted
256          * as a number.
257          * @memberof JXG.Math.Statistics
258          */
259         min: function (arr) {
260             return Math.min.apply(this, arr);
261         },
262 
263         /**
264          * Determines the lowest and the highest value from the given array.
265          * @param {Array} arr
266          * @returns {Array} The minimum value as the first and the maximum value as the second value.
267          * @memberof JXG.Math.Statistics
268          */
269         range: function (arr) {
270             return [this.min(arr), this.max(arr)];
271         },
272 
273         /**
274          * Determines the absolute value of every given value.
275          * @param {Array|Number} arr
276          * @returns {Array|Number}
277          * @memberof JXG.Math.Statistics
278          */
279         abs: function (arr) {
280             var i, len, res;
281 
282             if (Type.isArray(arr)) {
283                 if (arr.map) {
284                     res = arr.map(Math.abs);
285                 } else {
286                     len = arr.length;
287                     res = [];
288 
289                     for (i = 0; i < len; i++) {
290                         res[i] = Math.abs(arr[i]);
291                     }
292                 }
293             } else if (ArrayBuffer.isView(arr)) {
294                 res = arr.map(Math.abs);
295             } else {
296                 res = Math.abs(arr);
297             }
298             return res;
299         },
300 
301         /**
302          * Adds up two (sequences of) values. If one value is an array and the other one is a number the number
303          * is added to every element of the array. If two arrays are given and the lengths don't match the shortest
304          * length is taken.
305          * @param {Array|Number} arr1
306          * @param {Array|Number} arr2
307          * @returns {Array|Number}
308          * @memberof JXG.Math.Statistics
309          */
310         add: function (arr1, arr2) {
311             var i, len, res = [];
312 
313             arr1 = Type.evalSlider(arr1);
314             arr2 = Type.evalSlider(arr2);
315 
316             if (Type.isArray(arr1) && Type.isNumber(arr2)) {
317                 len = arr1.length;
318 
319                 for (i = 0; i < len; i++) {
320                     res[i] = arr1[i] + arr2;
321                 }
322             } else if (Type.isNumber(arr1) && Type.isArray(arr2)) {
323                 len = arr2.length;
324 
325                 for (i = 0; i < len; i++) {
326                     res[i] = arr1 + arr2[i];
327                 }
328             } else if (Type.isArray(arr1) && Type.isArray(arr2)) {
329                 len = Math.min(arr1.length, arr2.length);
330 
331                 for (i = 0; i < len; i++) {
332                     res[i] = arr1[i] + arr2[i];
333                 }
334             } else {
335                 res = arr1 + arr2;
336             }
337 
338             return res;
339         },
340 
341         /**
342          * Divides two (sequences of) values. If two arrays are given and the lengths don't match the shortest length
343          * is taken.
344          * @param {Array|Number} arr1 Dividend
345          * @param {Array|Number} arr2 Divisor
346          * @returns {Array|Number}
347          * @memberof JXG.Math.Statistics
348          */
349         div: function (arr1, arr2) {
350             var i, len, res = [];
351 
352             arr1 = Type.evalSlider(arr1);
353             arr2 = Type.evalSlider(arr2);
354 
355             if (Type.isArray(arr1) && Type.isNumber(arr2)) {
356                 len = arr1.length;
357 
358                 for (i = 0; i < len; i++) {
359                     res[i] = arr1[i] / arr2;
360                 }
361             } else if (Type.isNumber(arr1) && Type.isArray(arr2)) {
362                 len = arr2.length;
363 
364                 for (i = 0; i < len; i++) {
365                     res[i] = arr1 / arr2[i];
366                 }
367             } else if (Type.isArray(arr1) && Type.isArray(arr2)) {
368                 len = Math.min(arr1.length, arr2.length);
369 
370                 for (i = 0; i < len; i++) {
371                     res[i] = arr1[i] / arr2[i];
372                 }
373             } else {
374                 res = arr1 / arr2;
375             }
376 
377             return res;
378         },
379 
380         /**
381          * @function
382          * @deprecated Use {@link JXG.Math.Statistics.div} instead.
383          * @memberof JXG.Math.Statistics
384          */
385         divide: function () {
386             JXG.deprecated('Statistics.divide()', 'Statistics.div()');
387             Mat.Statistics.div.apply(Mat.Statistics, arguments);
388         },
389 
390         /**
391          * Divides two (sequences of) values and returns the remainder. If two arrays are given and the lengths don't
392          * match the shortest length is taken.
393          * @param {Array|Number} arr1 Dividend
394          * @param {Array|Number} arr2 Divisor
395          * @param {Boolean} [math=false] Mathematical mod or symmetric mod? Default is symmetric, the JavaScript <tt>%</tt> operator.
396          * @returns {Array|Number}
397          * @memberof JXG.Math.Statistics
398          */
399         mod: function (arr1, arr2, math) {
400             var i, len, res = [], mod = function (a, m) {
401                 return a % m;
402             };
403 
404             math = Type.def(math, false);
405 
406             if (math) {
407                 mod = Mat.mod;
408             }
409 
410             arr1 = Type.evalSlider(arr1);
411             arr2 = Type.evalSlider(arr2);
412 
413             if (Type.isArray(arr1) && Type.isNumber(arr2)) {
414                 len = arr1.length;
415 
416                 for (i = 0; i < len; i++) {
417                     res[i] = mod(arr1[i], arr2);
418                 }
419             } else if (Type.isNumber(arr1) && Type.isArray(arr2)) {
420                 len = arr2.length;
421 
422                 for (i = 0; i < len; i++) {
423                     res[i] = mod(arr1, arr2[i]);
424                 }
425             } else if (Type.isArray(arr1) && Type.isArray(arr2)) {
426                 len = Math.min(arr1.length, arr2.length);
427 
428                 for (i = 0; i < len; i++) {
429                     res[i] = mod(arr1[i], arr2[i]);
430                 }
431             } else {
432                 res = mod(arr1, arr2);
433             }
434 
435             return res;
436         },
437 
438         /**
439          * Multiplies two (sequences of) values. If one value is an array and the other one is a number the number
440          * is multiplied to every element of the array. If two arrays are given and the lengths don't match the shortest
441          * length is taken.
442          * @param {Array|Number} arr1
443          * @param {Array|Number} arr2
444          * @returns {Array|Number}
445          * @memberof JXG.Math.Statistics
446          */
447         multiply: function (arr1, arr2) {
448             var i, len, res = [];
449 
450             arr1 = Type.evalSlider(arr1);
451             arr2 = Type.evalSlider(arr2);
452 
453             if (Type.isArray(arr1) && Type.isNumber(arr2)) {
454                 len = arr1.length;
455 
456                 for (i = 0; i < len; i++) {
457                     res[i] = arr1[i] * arr2;
458                 }
459             } else if (Type.isNumber(arr1) && Type.isArray(arr2)) {
460                 len = arr2.length;
461 
462                 for (i = 0; i < len; i++) {
463                     res[i] = arr1 * arr2[i];
464                 }
465             } else if (Type.isArray(arr1) && Type.isArray(arr2)) {
466                 len = Math.min(arr1.length, arr2.length);
467 
468                 for (i = 0; i < len; i++) {
469                     res[i] = arr1[i] * arr2[i];
470                 }
471             } else {
472                 res = arr1 * arr2;
473             }
474 
475             return res;
476         },
477 
478         /**
479          * Subtracts two (sequences of) values. If two arrays are given and the lengths don't match the shortest
480          * length is taken.
481          * @param {Array|Number} arr1 Minuend
482          * @param {Array|Number} arr2 Subtrahend
483          * @returns {Array|Number}
484          * @memberof JXG.Math.Statistics
485          */
486         subtract: function (arr1, arr2) {
487             var i, len, res = [];
488 
489             arr1 = Type.evalSlider(arr1);
490             arr2 = Type.evalSlider(arr2);
491 
492             if (Type.isArray(arr1) && Type.isNumber(arr2)) {
493                 len = arr1.length;
494 
495                 for (i = 0; i < len; i++) {
496                     res[i] = arr1[i] - arr2;
497                 }
498             } else if (Type.isNumber(arr1) && Type.isArray(arr2)) {
499                 len = arr2.length;
500 
501                 for (i = 0; i < len; i++) {
502                     res[i] = arr1 - arr2[i];
503                 }
504             } else if (Type.isArray(arr1) && Type.isArray(arr2)) {
505                 len = Math.min(arr1.length, arr2.length);
506 
507                 for (i = 0; i < len; i++) {
508                     res[i] = arr1[i] - arr2[i];
509                 }
510             } else {
511                 res = arr1 - arr2;
512             }
513 
514             return res;
515         },
516 
517         /**
518          * The Theil-Sen estimator can be used to determine a more robust linear regression of a set of sample
519          * points than least squares regression in {@link JXG.Math.Numerics.regressionPolynomial}.
520          *
521          * If the function should be applied to an array a of points, a the coords array can be generated with
522          * JavaScript array.map:
523          *
524          * <pre>
525          * JXG.Math.Statistics.TheilSenRegression(a.map(el => el.coords));
526          * </pre>
527          *
528          * @param {Array} coords Array of {@link JXG.Coords}.
529          * @returns {Array} A stdform array of the regression line.
530          * @memberof JXG.Math.Statistics
531          *
532          * @example
533          * var board = JXG.JSXGraph.initBoard('jxgbox', { boundingbox: [-6,6,6,-6], axis : true });
534          * var a=[];
535          * a[0]=board.create('point', [0,0]);
536          * a[1]=board.create('point', [3,0]);
537          * a[2]=board.create('point', [0,3]);
538          *
539          * board.create('line', [
540          *     () => JXG.Math.Statistics.TheilSenRegression(a.map(el => el.coords))
541          *   ],
542          *   {strokeWidth:1, strokeColor:'black'});
543          *
544          * </pre><div id="JXG0a28be85-91c5-44d3-aae6-114e81217cf0" class="jxgbox" style="width: 300px; height: 300px;"></div>
545          * <script type="text/javascript">
546          *     (function() {
547          *         var board = JXG.JSXGraph.initBoard('JXG0a28be85-91c5-44d3-aae6-114e81217cf0',
548          *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
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          *     })();
561          *
562          * </script><pre>
563          *
564          */
565         TheilSenRegression: function (coords) {
566             var i, j,
567                 slopes = [],
568                 tmpslopes = [],
569                 yintercepts = [];
570 
571             for (i = 0; i < coords.length; i++) {
572                 tmpslopes.length = 0;
573 
574                 for (j = 0; j < coords.length; j++) {
575                     if (Math.abs(coords[j].usrCoords[1] - coords[i].usrCoords[1]) > Mat.eps) {
576                         tmpslopes[j] = (coords[j].usrCoords[2] - coords[i].usrCoords[2]) /
577                             (coords[j].usrCoords[1] - coords[i].usrCoords[1]);
578                     }
579                 }
580 
581                 slopes[i] = this.median(tmpslopes);
582                 yintercepts.push(coords[i].usrCoords[2] - slopes[i] * coords[i].usrCoords[1]);
583             }
584 
585             return [this.median(yintercepts), this.median(slopes), -1];
586         },
587 
588         /**
589          * Generate values of a standard normal random variable with the Marsaglia polar method, see
590          * https://en.wikipedia.org/wiki/Marsaglia_polar_method .
591          *
592          * @param {Number} mean mean value of the normal distribution
593          * @param {Number} stdDev standard deviation of the normal distribution
594          * @returns {Number} value of a standard normal random variable
595          */
596         generateGaussian: function (mean, stdDev) {
597             var u, v, s;
598 
599             if (this.hasSpare) {
600                 this.hasSpare = false;
601                 return this.spare * stdDev + mean;
602             }
603 
604             do {
605                 u = Math.random() * 2 - 1;
606                 v = Math.random() * 2 - 1;
607                 s = u * u + v * v;
608             } while (s >= 1 || s === 0);
609 
610             s = Math.sqrt(-2.0 * Math.log(s) / s);
611 
612             this.spare = v * s;
613             this.hasSpare = true;
614             return mean + stdDev * u * s;
615         }
616     };
617 
618     return Mat.Statistics;
619 });
620