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