1 /* 2 Copyright 2008-2024 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 /** 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 * {@link https://en.wikipedia.org/wiki/Marsaglia_polar_method}. 607 * See also D. E. Knuth, The art of computer programming, vol 2, p. 117. 608 * 609 * @param {Number} mean mean value of the normal distribution 610 * @param {Number} stdDev standard deviation of the normal distribution 611 * @returns {Number} value of a standard normal random variable 612 * @memberof JXG.Math.Statistics 613 */ 614 generateGaussian: function (mean, stdDev) { 615 var u, v, s; 616 617 if (this.hasSpare) { 618 this.hasSpare = false; 619 return this.spare * stdDev + mean; 620 } 621 622 do { 623 u = Math.random() * 2 - 1; 624 v = Math.random() * 2 - 1; 625 s = u * u + v * v; 626 } while (s >= 1 || s === 0); 627 628 s = Math.sqrt((-2.0 * Math.log(s)) / s); 629 630 this.spare = v * s; 631 this.hasSpare = true; 632 return mean + stdDev * u * s; 633 }, 634 635 /** 636 * Generate value of a standard normal random variable with given mean and standard deviation. 637 * Alias for {@link JXG.Math.Statistics#generateGaussian} 638 * 639 * @param {Number} mean 640 * @param {Number} stdDev 641 * @returns Number 642 * @memberof JXG.Math.Statistics 643 * @see JXG.Math.Statistics#generateGaussian 644 */ 645 randomNormal: function (mean, stdDev) { 646 return this.generateGaussian(mean, stdDev); 647 }, 648 649 /** 650 * Generate value of a uniform distributed random variable in the interval [a, b]. 651 * @param {Number} a 652 * @param {Number} b 653 * @returns Number 654 * @memberof JXG.Math.Statistics 655 */ 656 randomUniform: function (a, b) { 657 return Math.random() * (b - a) + a; 658 }, 659 660 /** 661 * Generate value of a random variable with exponential distribution, i.e. 662 * <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>. 663 * See {@link https://en.wikipedia.org/wiki/Exponential_distribution}. 664 * Algorithm: D.E. Knuth, TAOCP 2, p. 128. 665 * 666 * @param {Number} lambda <i>> 0</i> 667 * @returns Number 668 * @memberof JXG.Math.Statistics 669 */ 670 randomExponential: function (lbda) { 671 var u; 672 673 // Knuth, TAOCP 2, p 128 674 // See https://en.wikipedia.org/wiki/Exponential_distribution 675 if (lbda <= 0) { 676 return NaN; 677 } 678 679 do { 680 u = Math.random(); 681 } while (u === 0); 682 683 return -Math.log(u) / lbda; 684 }, 685 686 /** 687 * Generate value of a random variable with gamma distribution of order alpha. 688 * See {@link https://en.wikipedia.org/wiki/Gamma_distribution}. 689 * Algorithm: D.E. Knuth, TAOCP 2, p. 129. 690 691 * @param {Number} a shape, <i> > 0</i> 692 * @param {Number} [b=1] scale, <i> > 0</i> 693 * @param {Number} [t=0] threshold 694 * @returns Number 695 * @memberof JXG.Math.Statistics 696 */ 697 randomGamma: function (a, b, t) { 698 var u, v, x, y, 699 p, q; 700 701 if (a <= 0) { 702 return NaN; 703 } 704 705 b = b || 1; 706 t = t || 0; 707 708 if (a === 1) { 709 return b * this.randomExponential(1) + t; 710 } 711 712 if (a < 1) { 713 // Method by Ahrens 714 // Knuth, TAOCP 2, Ex. 16, p 551 715 p = Math.E / (a + Math.E); 716 717 do { 718 u = Math.random(); 719 do { 720 v = Math.random(); 721 } while (v === 0); 722 if (u < p) { 723 x = Math.pow(v, 1 / a); 724 q = Math.exp(-x); 725 } else { 726 x = 1 - Math.log(v); 727 q = Math.pow(x, a - 1); 728 } 729 u = Math.random(); 730 } while (u >= q); 731 return b * x + t; 732 } 733 734 // a > 1 735 // Knuth, TAOCP 2, p 129 736 do { 737 y = Math.tan(Math.PI * Math.random()); 738 x = Math.sqrt(2 * a - 1) * y + a - 1; 739 if (x > 0) { 740 v = Math.random(); 741 } else { 742 continue; 743 } 744 } while (x <= 0.0 || v > (1 + y * y) * Math.exp( (a - 1) * Math.log(x / (a-1)) - Math.sqrt(2 * a - 1) * y)); 745 746 return b * x + t; 747 }, 748 749 /** 750 * Generate value of a random variable with beta distribution with shape parameters alpha and beta. 751 * See {@link https://en.wikipedia.org/wiki/Beta_distribution}. 752 * 753 * @param {Number} alpha <i>> 0</i> 754 * @param {Number} beta <i>> 0</i> 755 * @returns Number 756 * @memberof JXG.Math.Statistics 757 */ 758 randomBeta: function (a, b) { 759 // Knuth, TAOCP 2, p 129 760 var x1, x2, x; 761 762 if (a <= 0 || b <= 0) { 763 return NaN; 764 } 765 766 x1 = this.randomGamma(a); 767 x2 = this.randomGamma(b); 768 x = x1 / (x1 + x2); 769 return x; 770 }, 771 772 /** 773 * Generate value of a random variable with chi-square distribution with k degrees of freedom. 774 * See {@link https://en.wikipedia.org/wiki/Chi-squared_distribution}. 775 * 776 * @param {Number} k <i>> 0</i> 777 * @returns Number 778 * @memberof JXG.Math.Statistics 779 */ 780 randomChisquare: function (nu) { 781 // Knuth, TAOCP 2, p 130 782 783 if (nu <= 0) { 784 return NaN; 785 } 786 787 return 2 * this.randomGamma(nu * 0.5); 788 }, 789 790 /** 791 * Generate value of a random variable with F-distribution with d<sub>1</sub> and d<sub>2</sub> degrees of freedom. 792 * See {@link https://en.wikipedia.org/wiki/F-distribution}. 793 * @param {Number} d1 <i>> 0</i> 794 * @param {Number} d2 <i>> 0</i> 795 * @returns Number 796 * @memberof JXG.Math.Statistics 797 */ 798 randomF: function (nu1, nu2) { 799 // Knuth, TAOCP 2, p 130 800 var y1, y2; 801 802 if (nu1 <= 0 || nu2 <= 0) { 803 return NaN; 804 } 805 806 y1 = this.randomChisquare(nu1); 807 y2 = this.randomChisquare(nu2); 808 809 return (y1 * nu2) / (y2 * nu1); 810 }, 811 812 /** 813 * Generate value of a random variable with Students-t-distribution with ν degrees of freedom. 814 * See {@link https://en.wikipedia.org/wiki/Student%27s_t-distribution}. 815 * @param {Number} nu <i>> 0</i> 816 * @returns Number 817 * @memberof JXG.Math.Statistics 818 */ 819 randomT: function (nu) { 820 // Knuth, TAOCP 2, p 130 821 var y1, y2; 822 823 if (nu <= 0) { 824 return NaN; 825 } 826 827 y1 = this.randomNormal(0, 1); 828 y2 = this.randomChisquare(nu); 829 830 return y1 / Math.sqrt(y2 / nu); 831 }, 832 833 /** 834 * Generate values for a random variable in binomial distribution with parameters <i>n</i> and <i>p</i>. 835 * See {@link https://en.wikipedia.org/wiki/Binomial_distribution}. 836 * It uses algorithm BG from {@link https://dl.acm.org/doi/pdf/10.1145/42372.42381}. 837 * 838 * @param {Number} n Number of trials (n >= 0) 839 * @param {Number} p Propability (0 <= p <= 1) 840 * @returns Number Integer value of a random variable in binomial distribution 841 * @memberof JXG.Math.Statistics 842 * 843 * @example 844 * console.log(JXG.Mat.Statistics.generateBinomial(100,0.1)); 845 * // Possible output: 18 846 * 847 */ 848 randomBinomial: function (n, p) { 849 var x, y, c, 850 a, b, N1; 851 852 if (p < 0 || p > 1 || n < 0) { 853 return NaN; 854 } 855 856 // Edge cases 857 if (p === 0) { 858 return 0; 859 } 860 if (p === 1) { 861 return n; 862 } 863 864 // Now, we can assume 0 < p < 1. 865 866 // Fast path for common cases 867 if (n === 0) { 868 return 0; 869 } 870 if (n === 1) { 871 return ((Math.random() < p) ? 1 : 0); 872 } 873 874 // Exploit symmetry 875 if (p > 0.5) { 876 return n - this.randomBinomial(n, 1 - p); 877 } 878 879 // General case: n > 1, p <= 0.5 880 if (n < 100) { 881 // n small: 882 // Algorithm BG (Devroye) from: 883 // https://dl.acm.org/doi/pdf/10.1145/42372.42381 884 // Time O(np) so suitable for np small only. 885 x = -1; 886 y = 0; 887 888 c = Math.log(1 - p); 889 if (c === 0) { 890 return 0; 891 } 892 893 do { 894 x += 1; 895 y += Math.floor(Math.log(Math.random()) / c) + 1; 896 } while (y < n); 897 } else { 898 // n large: 899 // Knuth, TAOCP 2, p 131 900 a = 1 + Math.floor(n * 0.5); 901 b = n - a + 1; 902 x = this.randomBeta(a, b); 903 if (x >= p) { 904 N1 = this.randomBinomial(a - 1, p / x); 905 x = N1; 906 } else { 907 N1 = this.randomBinomial(b - 1, (p - x) / (1 - x)); 908 x = a + N1; 909 } 910 } 911 return x; 912 }, 913 914 /** 915 * Generate values for a random variable in geometric distribution with propability <i>p</i>. 916 * See {@link https://en.wikipedia.org/wiki/Geometric_distribution}. 917 * 918 * @param {Number} p (0 <= p <= 1) 919 * @returns Number 920 * @memberof JXG.Math.Statistics 921 */ 922 randomGeometric: function(p) { 923 var u; 924 925 if (p < 0 || p > 1) { 926 return NaN; 927 } 928 // Knuth, TAOCP 2, p 131 929 u = Math.random(); 930 931 return Math.ceil(Math.log(u) / Math.log(1 - p)); 932 }, 933 934 /** 935 * Generate values for a random variable in Poisson distribution with mean <i>mu</i>. 936 * See {@link https://en.wikipedia.org/wiki/Poisson_distribution}. 937 * 938 * @param {Number} mu (0 < mu) 939 * @returns Number 940 * @memberof JXG.Math.Statistics 941 */ 942 randomPoisson: function (mu) { 943 var e = Math.exp(-mu), 944 N, 945 m = 0, 946 u = 1, 947 x, 948 alpha = 7 / 8; 949 950 if (mu <= 0) { 951 return NaN; 952 } 953 954 // Knuth, TAOCP 2, p 132 955 if (mu < 10) { 956 do { 957 u *= Math.random(); 958 m += 1; 959 } while (u > e); 960 N = m - 1; 961 } else { 962 m = Math.floor(alpha * mu); 963 x = this.randomGamma(m); 964 if (x < mu) { 965 N = m + this.randomPoisson(mu - x); 966 } else { 967 N = this.randomBinomial(m - 1, mu / x); 968 } 969 } 970 return N; 971 }, 972 973 /** 974 * Generate values for a random variable in hypergeometric distribution. 975 * Samples are drawn from a hypergeometric distribution with specified parameters, <i>good</i> (ways to make a good selection), 976 * <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>). 977 * <p> 978 * Naive implementation with runtime <i>O(samples)</i>. 979 * 980 * @param {Number} good ways to make a good selection 981 * @param {Number} bad ways to make a bad selection 982 * @param {Number} samples number of items sampled 983 * @returns 984 * @memberof JXG.Math.Statistics 985 */ 986 randomHypergeometric: function (good, bad, k) { 987 var i, u, 988 x = 0, 989 // kk, 990 // n = good + bad, 991 d1 = good + bad - k, 992 d2 = Math.min(good, bad), 993 y = d2; 994 995 if (good < 1 || bad < 1 || k > good + bad) { 996 return NaN; 997 } 998 999 // Naive method 1000 // kk = Math.min(k, n - k); 1001 // for (i = 0; i < k; i ++) { 1002 // u = Math.random(); 1003 // if (n * u <= good) { 1004 // x += 1; 1005 // if (x === good) { 1006 // return x; 1007 // } 1008 // good -= 1; 1009 // } 1010 // n -= 1; 1011 // } 1012 // return x; 1013 1014 // Implementation from 1015 // Monte Carlo by George S. Fishman 1016 // https://link.springer.com/book/10.1007/978-1-4757-2553-7 1017 // page 218 1018 // 1019 i = k; 1020 while (y * i > 0) { 1021 u = Math.random(); 1022 y -= Math.floor(u + y / (d1 + i)); 1023 i -= 1; 1024 } 1025 x = d2 - y; 1026 if (good <= bad) { 1027 return x; 1028 } else { 1029 return k - x; 1030 } 1031 }, 1032 1033 /** 1034 * Compute the histogram of a dataset. 1035 * Optional parameters can be supplied through a JavaScript object 1036 * with the following default values: 1037 * <pre> 1038 * { 1039 * bins: 10, // Number of bins 1040 * range: false, // false or array. The lower and upper range of the bins. 1041 * // If not provided, range is simply [min(x), max(x)]. 1042 * // Values outside the range are ignored. 1043 * density: false, // If true, normalize the counts by dividing by sum(counts) 1044 * cumulative: false 1045 * } 1046 * </pre> 1047 * The function returns an array containing two arrays. The first array is of length bins+1 1048 * containing the start values of the bins. The last entry contains the end values of the last bin. 1049 * <p> 1050 * The second array contains the counts of each bin. 1051 * @param {Array} x 1052 * @param {Object} opt Optional parameters 1053 * @returns Array [bin, counts] Array bins contains start values of bins, array counts contains 1054 * the number of entries of x which are contained in each bin. 1055 * @memberof JXG.Math.Statistics 1056 * 1057 * @example 1058 * var curve = board.create('curve', [[], []]); 1059 * curve.updateDataArray = function () { 1060 * var i, res, x = []; 1061 * 1062 * for (i = 0; i < 5000; i++) { 1063 * x.push(JXG.Math.Statistics.randomGamma(2)); 1064 * } 1065 * res = JXG.Math.Statistics.histogram(x, { bins: 50, density: true, cumulative: false, range: [-5, 5] }); 1066 * this.dataX = res[1]; 1067 * this.dataY = res[0]; 1068 * }; 1069 * board.update(); 1070 * 1071 * </pre><div id="JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302" class="jxgbox" style="width: 300px; height: 300px;"></div> 1072 * <script type="text/javascript"> 1073 * (function() { 1074 * var board = JXG.JSXGraph.initBoard('JXGda56df4d-a5a5-4c87-9ffc-9bbc1b512302', 1075 * {boundingbox: [-1, 3, 6, -1], axis: true, showcopyright: false, shownavigation: false}); 1076 * var curve = board.create('curve', [[], []]); 1077 * curve.updateDataArray = function () { 1078 * var i, res, x = []; 1079 * 1080 * for (i = 0; i < 5000; i++) { 1081 * x.push(JXG.Math.Statistics.randomGamma(2)); 1082 * } 1083 * res = JXG.Math.Statistics.histogram(x, { bins: 50, density: true, cumulative: false, range: [-5, 5] }); 1084 * this.dataX = res[1]; 1085 * this.dataY = res[0]; 1086 * }; 1087 * board.update(); 1088 * })(); 1089 * 1090 * </script><pre> 1091 * 1092 */ 1093 histogram: function (x, opt) { 1094 var i, le, k, 1095 mi, ma, num_bins, delta, 1096 range, 1097 s, 1098 counts = [], 1099 bins = []; 1100 1101 // Evaluate number of bins 1102 num_bins = opt.bins || 10; 1103 1104 // Evaluate range 1105 range = opt.range || false; 1106 if (range === false) { 1107 mi = Math.min.apply(null, x); 1108 ma = Math.max.apply(null, x); 1109 } else { 1110 mi = range[0]; 1111 ma = range[1]; 1112 } 1113 1114 // Set uniform delta 1115 if (num_bins > 0) { 1116 delta = (ma - mi) / (num_bins - 1); 1117 } else { 1118 delta = 0; 1119 } 1120 1121 // Set the bins and init the counts array 1122 for (i = 0; i < num_bins; i++) { 1123 counts.push(0); 1124 bins.push(mi + i * delta); 1125 } 1126 // bins.push(ma); 1127 1128 // Determine the counts 1129 le = x.length; 1130 for (i = 0; i < le; i++) { 1131 k = Math.floor((x[i] - mi) / delta); 1132 if (k >= 0 && k < num_bins) { 1133 counts[k] += 1; 1134 } 1135 } 1136 1137 // Normalize if density===true 1138 if (opt.density) { 1139 s = JXG.Math.Statistics.sum(counts); 1140 for (i = 0; i < num_bins; i++) { 1141 // counts[i] /= (s * delta); 1142 counts[i] /= s; 1143 } 1144 } 1145 1146 // Cumulative counts 1147 if (opt.cumulative) { 1148 for (i = 1; i < num_bins; i++) { 1149 counts[i] += counts[i - 1]; 1150 } 1151 } 1152 1153 return [counts, bins]; 1154 } 1155 }; 1156 1157 export default Mat.Statistics; 1158