      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