1 /* 2 Copyright 2008-2021 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