1 | (function (window, factory) {
|
2 | if (typeof exports === 'object') {
|
3 | module.exports = factory();
|
4 | } else if (typeof define === 'function' && define.amd) {
|
5 | define(factory);
|
6 | } else {
|
7 | window.jStat = factory();
|
8 | }
|
9 | })(this, function () {
|
10 | var jStat = (function(Math, undefined) {
|
11 |
|
12 | console.warn('The npm package jStat is no longer maintained. Instead use jstat (lowercase).');
|
13 | console.warn('Visit https://www.npmjs.com/package/jstat for more information.');
|
14 |
|
15 |
|
16 | var concat = Array.prototype.concat;
|
17 | var slice = Array.prototype.slice;
|
18 | var toString = Object.prototype.toString;
|
19 |
|
20 |
|
21 |
|
22 | function calcRdx(n, m) {
|
23 | var val = n > m ? n : m;
|
24 | return Math.pow(10,
|
25 | 17 - ~~(Math.log(((val > 0) ? val : -val)) * Math.LOG10E));
|
26 | }
|
27 |
|
28 |
|
29 | var isArray = Array.isArray || function isArray(arg) {
|
30 | return toString.call(arg) === '[object Array]';
|
31 | };
|
32 |
|
33 |
|
34 | function isFunction(arg) {
|
35 | return toString.call(arg) === '[object Function]';
|
36 | }
|
37 |
|
38 |
|
39 | function isNumber(arg) {
|
40 | return typeof arg === 'number' && arg === arg;
|
41 | }
|
42 |
|
43 |
|
44 |
|
45 | function toVector(arr) {
|
46 | return concat.apply([], arr);
|
47 | }
|
48 |
|
49 |
|
50 |
|
51 | function jStat() {
|
52 | return new jStat._init(arguments);
|
53 | }
|
54 |
|
55 |
|
56 |
|
57 | jStat.fn = jStat.prototype;
|
58 |
|
59 |
|
60 |
|
61 |
|
62 | jStat._init = function _init(args) {
|
63 |
|
64 | if (isArray(args[0])) {
|
65 |
|
66 | if (isArray(args[0][0])) {
|
67 |
|
68 | if (isFunction(args[1]))
|
69 | args[0] = jStat.map(args[0], args[1]);
|
70 |
|
71 | for (var i = 0; i < args[0].length; i++)
|
72 | this[i] = args[0][i];
|
73 | this.length = args[0].length;
|
74 |
|
75 |
|
76 | } else {
|
77 | this[0] = isFunction(args[1]) ? jStat.map(args[0], args[1]) : args[0];
|
78 | this.length = 1;
|
79 | }
|
80 |
|
81 |
|
82 | } else if (isNumber(args[0])) {
|
83 | this[0] = jStat.seq.apply(null, args);
|
84 | this.length = 1;
|
85 |
|
86 |
|
87 | } else if (args[0] instanceof jStat) {
|
88 |
|
89 | return jStat(args[0].toArray());
|
90 |
|
91 |
|
92 |
|
93 |
|
94 | } else {
|
95 | this[0] = [];
|
96 | this.length = 1;
|
97 | }
|
98 |
|
99 | return this;
|
100 | };
|
101 | jStat._init.prototype = jStat.prototype;
|
102 | jStat._init.constructor = jStat;
|
103 |
|
104 |
|
105 |
|
106 |
|
107 | jStat.utils = {
|
108 | calcRdx: calcRdx,
|
109 | isArray: isArray,
|
110 | isFunction: isFunction,
|
111 | isNumber: isNumber,
|
112 | toVector: toVector
|
113 | };
|
114 |
|
115 |
|
116 | jStat._random_fn = Math.random;
|
117 | jStat.setRandom = function setRandom(fn) {
|
118 | if (typeof fn !== 'function')
|
119 | throw new TypeError('fn is not a function');
|
120 | jStat._random_fn = fn;
|
121 | };
|
122 |
|
123 |
|
124 |
|
125 |
|
126 | jStat.extend = function extend(obj) {
|
127 | var i, j;
|
128 |
|
129 | if (arguments.length === 1) {
|
130 | for (j in obj)
|
131 | jStat[j] = obj[j];
|
132 | return this;
|
133 | }
|
134 |
|
135 | for (i = 1; i < arguments.length; i++) {
|
136 | for (j in arguments[i])
|
137 | obj[j] = arguments[i][j];
|
138 | }
|
139 |
|
140 | return obj;
|
141 | };
|
142 |
|
143 |
|
144 |
|
145 | jStat.rows = function rows(arr) {
|
146 | return arr.length || 1;
|
147 | };
|
148 |
|
149 |
|
150 |
|
151 | jStat.cols = function cols(arr) {
|
152 | return arr[0].length || 1;
|
153 | };
|
154 |
|
155 |
|
156 |
|
157 | jStat.dimensions = function dimensions(arr) {
|
158 | return {
|
159 | rows: jStat.rows(arr),
|
160 | cols: jStat.cols(arr)
|
161 | };
|
162 | };
|
163 |
|
164 |
|
165 |
|
166 | jStat.row = function row(arr, index) {
|
167 | if (isArray(index)) {
|
168 | return index.map(function(i) {
|
169 | return jStat.row(arr, i);
|
170 | })
|
171 | }
|
172 | return arr[index];
|
173 | };
|
174 |
|
175 |
|
176 |
|
177 |
|
178 | jStat.rowa = function rowa(arr, i) {
|
179 | return jStat.row(arr, i);
|
180 | };
|
181 |
|
182 |
|
183 |
|
184 |
|
185 | jStat.col = function col(arr, index) {
|
186 | if (isArray(index)) {
|
187 | var submat = jStat.arange(arr.length).map(function() {
|
188 | return new Array(index.length);
|
189 | });
|
190 | index.forEach(function(ind, i){
|
191 | jStat.arange(arr.length).forEach(function(j) {
|
192 | submat[j][i] = arr[j][ind];
|
193 | });
|
194 | });
|
195 | return submat;
|
196 | }
|
197 | var column = new Array(arr.length);
|
198 | for (var i = 0; i < arr.length; i++)
|
199 | column[i] = [arr[i][index]];
|
200 | return column;
|
201 | };
|
202 |
|
203 |
|
204 |
|
205 |
|
206 | jStat.cola = function cola(arr, i) {
|
207 | return jStat.col(arr, i).map(function(a){ return a[0] });
|
208 | };
|
209 |
|
210 |
|
211 |
|
212 | jStat.diag = function diag(arr) {
|
213 | var nrow = jStat.rows(arr);
|
214 | var res = new Array(nrow);
|
215 | for (var row = 0; row < nrow; row++)
|
216 | res[row] = [arr[row][row]];
|
217 | return res;
|
218 | };
|
219 |
|
220 |
|
221 |
|
222 | jStat.antidiag = function antidiag(arr) {
|
223 | var nrow = jStat.rows(arr) - 1;
|
224 | var res = new Array(nrow);
|
225 | for (var i = 0; nrow >= 0; nrow--, i++)
|
226 | res[i] = [arr[i][nrow]];
|
227 | return res;
|
228 | };
|
229 |
|
230 |
|
231 | jStat.transpose = function transpose(arr) {
|
232 | var obj = [];
|
233 | var objArr, rows, cols, j, i;
|
234 |
|
235 |
|
236 | if (!isArray(arr[0]))
|
237 | arr = [arr];
|
238 |
|
239 | rows = arr.length;
|
240 | cols = arr[0].length;
|
241 |
|
242 | for (i = 0; i < cols; i++) {
|
243 | objArr = new Array(rows);
|
244 | for (j = 0; j < rows; j++)
|
245 | objArr[j] = arr[j][i];
|
246 | obj.push(objArr);
|
247 | }
|
248 |
|
249 |
|
250 | return obj.length === 1 ? obj[0] : obj;
|
251 | };
|
252 |
|
253 |
|
254 |
|
255 |
|
256 | jStat.map = function map(arr, func, toAlter) {
|
257 | var row, nrow, ncol, res, col;
|
258 |
|
259 | if (!isArray(arr[0]))
|
260 | arr = [arr];
|
261 |
|
262 | nrow = arr.length;
|
263 | ncol = arr[0].length;
|
264 | res = toAlter ? arr : new Array(nrow);
|
265 |
|
266 | for (row = 0; row < nrow; row++) {
|
267 |
|
268 | if (!res[row])
|
269 | res[row] = new Array(ncol);
|
270 | for (col = 0; col < ncol; col++)
|
271 | res[row][col] = func(arr[row][col], row, col);
|
272 | }
|
273 |
|
274 | return res.length === 1 ? res[0] : res;
|
275 | };
|
276 |
|
277 |
|
278 |
|
279 | jStat.cumreduce = function cumreduce(arr, func, toAlter) {
|
280 | var row, nrow, ncol, res, col;
|
281 |
|
282 | if (!isArray(arr[0]))
|
283 | arr = [arr];
|
284 |
|
285 | nrow = arr.length;
|
286 | ncol = arr[0].length;
|
287 | res = toAlter ? arr : new Array(nrow);
|
288 |
|
289 | for (row = 0; row < nrow; row++) {
|
290 |
|
291 | if (!res[row])
|
292 | res[row] = new Array(ncol);
|
293 | if (ncol > 0)
|
294 | res[row][0] = arr[row][0];
|
295 | for (col = 1; col < ncol; col++)
|
296 | res[row][col] = func(res[row][col-1], arr[row][col]);
|
297 | }
|
298 | return res.length === 1 ? res[0] : res;
|
299 | };
|
300 |
|
301 |
|
302 |
|
303 | jStat.alter = function alter(arr, func) {
|
304 | return jStat.map(arr, func, true);
|
305 | };
|
306 |
|
307 |
|
308 |
|
309 | jStat.create = function create(rows, cols, func) {
|
310 | var res = new Array(rows);
|
311 | var i, j;
|
312 |
|
313 | if (isFunction(cols)) {
|
314 | func = cols;
|
315 | cols = rows;
|
316 | }
|
317 |
|
318 | for (i = 0; i < rows; i++) {
|
319 | res[i] = new Array(cols);
|
320 | for (j = 0; j < cols; j++)
|
321 | res[i][j] = func(i, j);
|
322 | }
|
323 |
|
324 | return res;
|
325 | };
|
326 |
|
327 |
|
328 | function retZero() { return 0; }
|
329 |
|
330 |
|
331 |
|
332 | jStat.zeros = function zeros(rows, cols) {
|
333 | if (!isNumber(cols))
|
334 | cols = rows;
|
335 | return jStat.create(rows, cols, retZero);
|
336 | };
|
337 |
|
338 |
|
339 | function retOne() { return 1; }
|
340 |
|
341 |
|
342 |
|
343 | jStat.ones = function ones(rows, cols) {
|
344 | if (!isNumber(cols))
|
345 | cols = rows;
|
346 | return jStat.create(rows, cols, retOne);
|
347 | };
|
348 |
|
349 |
|
350 |
|
351 | jStat.rand = function rand(rows, cols) {
|
352 | if (!isNumber(cols))
|
353 | cols = rows;
|
354 | return jStat.create(rows, cols, jStat._random_fn);
|
355 | };
|
356 |
|
357 |
|
358 | function retIdent(i, j) { return i === j ? 1 : 0; }
|
359 |
|
360 |
|
361 |
|
362 | jStat.identity = function identity(rows, cols) {
|
363 | if (!isNumber(cols))
|
364 | cols = rows;
|
365 | return jStat.create(rows, cols, retIdent);
|
366 | };
|
367 |
|
368 |
|
369 |
|
370 | jStat.symmetric = function symmetric(arr) {
|
371 | var size = arr.length;
|
372 | var row, col;
|
373 |
|
374 | if (arr.length !== arr[0].length)
|
375 | return false;
|
376 |
|
377 | for (row = 0; row < size; row++) {
|
378 | for (col = 0; col < size; col++)
|
379 | if (arr[col][row] !== arr[row][col])
|
380 | return false;
|
381 | }
|
382 |
|
383 | return true;
|
384 | };
|
385 |
|
386 |
|
387 |
|
388 | jStat.clear = function clear(arr) {
|
389 | return jStat.alter(arr, retZero);
|
390 | };
|
391 |
|
392 |
|
393 |
|
394 | jStat.seq = function seq(min, max, length, func) {
|
395 | if (!isFunction(func))
|
396 | func = false;
|
397 |
|
398 | var arr = [];
|
399 | var hival = calcRdx(min, max);
|
400 | var step = (max * hival - min * hival) / ((length - 1) * hival);
|
401 | var current = min;
|
402 | var cnt;
|
403 |
|
404 |
|
405 |
|
406 | for (cnt = 0;
|
407 | current <= max && cnt < length;
|
408 | cnt++, current = (min * hival + step * hival * cnt) / hival) {
|
409 | arr.push((func ? func(current, cnt) : current));
|
410 | }
|
411 |
|
412 | return arr;
|
413 | };
|
414 |
|
415 |
|
416 |
|
417 |
|
418 |
|
419 | jStat.arange = function arange(start, end, step) {
|
420 | var rl = [];
|
421 | var i;
|
422 | step = step || 1;
|
423 | if (end === undefined) {
|
424 | end = start;
|
425 | start = 0;
|
426 | }
|
427 | if (start === end || step === 0) {
|
428 | return [];
|
429 | }
|
430 | if (start < end && step < 0) {
|
431 | return [];
|
432 | }
|
433 | if (start > end && step > 0) {
|
434 | return [];
|
435 | }
|
436 | if (step > 0) {
|
437 | for (i = start; i < end; i += step) {
|
438 | rl.push(i);
|
439 | }
|
440 | } else {
|
441 | for (i = start; i > end; i += step) {
|
442 | rl.push(i);
|
443 | }
|
444 | }
|
445 | return rl;
|
446 | };
|
447 |
|
448 |
|
449 |
|
450 |
|
451 |
|
452 |
|
453 | jStat.slice = (function(){
|
454 | function _slice(list, start, end, step) {
|
455 |
|
456 | var i;
|
457 | var rl = [];
|
458 | var length = list.length;
|
459 | if (start === undefined && end === undefined && step === undefined) {
|
460 | return jStat.copy(list);
|
461 | }
|
462 |
|
463 | start = start || 0;
|
464 | end = end || list.length;
|
465 | start = start >= 0 ? start : length + start;
|
466 | end = end >= 0 ? end : length + end;
|
467 | step = step || 1;
|
468 | if (start === end || step === 0) {
|
469 | return [];
|
470 | }
|
471 | if (start < end && step < 0) {
|
472 | return [];
|
473 | }
|
474 | if (start > end && step > 0) {
|
475 | return [];
|
476 | }
|
477 | if (step > 0) {
|
478 | for (i = start; i < end; i += step) {
|
479 | rl.push(list[i]);
|
480 | }
|
481 | } else {
|
482 | for (i = start; i > end;i += step) {
|
483 | rl.push(list[i]);
|
484 | }
|
485 | }
|
486 | return rl;
|
487 | }
|
488 |
|
489 | function slice(list, rcSlice) {
|
490 | var colSlice, rowSlice;
|
491 | rcSlice = rcSlice || {};
|
492 | if (isNumber(rcSlice.row)) {
|
493 | if (isNumber(rcSlice.col))
|
494 | return list[rcSlice.row][rcSlice.col];
|
495 | var row = jStat.rowa(list, rcSlice.row);
|
496 | colSlice = rcSlice.col || {};
|
497 | return _slice(row, colSlice.start, colSlice.end, colSlice.step);
|
498 | }
|
499 |
|
500 | if (isNumber(rcSlice.col)) {
|
501 | var col = jStat.cola(list, rcSlice.col);
|
502 | rowSlice = rcSlice.row || {};
|
503 | return _slice(col, rowSlice.start, rowSlice.end, rowSlice.step);
|
504 | }
|
505 |
|
506 | rowSlice = rcSlice.row || {};
|
507 | colSlice = rcSlice.col || {};
|
508 | var rows = _slice(list, rowSlice.start, rowSlice.end, rowSlice.step);
|
509 | return rows.map(function(row) {
|
510 | return _slice(row, colSlice.start, colSlice.end, colSlice.step);
|
511 | });
|
512 | }
|
513 |
|
514 | return slice;
|
515 | }());
|
516 |
|
517 |
|
518 |
|
519 |
|
520 |
|
521 | jStat.sliceAssign = function sliceAssign(A, rcSlice, B) {
|
522 | var nl, ml;
|
523 | if (isNumber(rcSlice.row)) {
|
524 | if (isNumber(rcSlice.col))
|
525 | return A[rcSlice.row][rcSlice.col] = B;
|
526 | rcSlice.col = rcSlice.col || {};
|
527 | rcSlice.col.start = rcSlice.col.start || 0;
|
528 | rcSlice.col.end = rcSlice.col.end || A[0].length;
|
529 | rcSlice.col.step = rcSlice.col.step || 1;
|
530 | nl = jStat.arange(rcSlice.col.start,
|
531 | Math.min(A.length, rcSlice.col.end),
|
532 | rcSlice.col.step);
|
533 | var m = rcSlice.row;
|
534 | nl.forEach(function(n, i) {
|
535 | A[m][n] = B[i];
|
536 | });
|
537 | return A;
|
538 | }
|
539 |
|
540 | if (isNumber(rcSlice.col)) {
|
541 | rcSlice.row = rcSlice.row || {};
|
542 | rcSlice.row.start = rcSlice.row.start || 0;
|
543 | rcSlice.row.end = rcSlice.row.end || A.length;
|
544 | rcSlice.row.step = rcSlice.row.step || 1;
|
545 | ml = jStat.arange(rcSlice.row.start,
|
546 | Math.min(A[0].length, rcSlice.row.end),
|
547 | rcSlice.row.step);
|
548 | var n = rcSlice.col;
|
549 | ml.forEach(function(m, j) {
|
550 | A[m][n] = B[j];
|
551 | });
|
552 | return A;
|
553 | }
|
554 |
|
555 | if (B[0].length === undefined) {
|
556 | B = [B];
|
557 | }
|
558 | rcSlice.row.start = rcSlice.row.start || 0;
|
559 | rcSlice.row.end = rcSlice.row.end || A.length;
|
560 | rcSlice.row.step = rcSlice.row.step || 1;
|
561 | rcSlice.col.start = rcSlice.col.start || 0;
|
562 | rcSlice.col.end = rcSlice.col.end || A[0].length;
|
563 | rcSlice.col.step = rcSlice.col.step || 1;
|
564 | ml = jStat.arange(rcSlice.row.start,
|
565 | Math.min(A.length, rcSlice.row.end),
|
566 | rcSlice.row.step);
|
567 | nl = jStat.arange(rcSlice.col.start,
|
568 | Math.min(A[0].length, rcSlice.col.end),
|
569 | rcSlice.col.step);
|
570 | ml.forEach(function(m, i) {
|
571 | nl.forEach(function(n, j) {
|
572 | A[m][n] = B[i][j];
|
573 | });
|
574 | });
|
575 | return A;
|
576 | };
|
577 |
|
578 |
|
579 |
|
580 |
|
581 | jStat.diagonal = function diagonal(diagArray) {
|
582 | var mat = jStat.zeros(diagArray.length, diagArray.length);
|
583 | diagArray.forEach(function(t, i) {
|
584 | mat[i][i] = t;
|
585 | });
|
586 | return mat;
|
587 | };
|
588 |
|
589 |
|
590 |
|
591 | jStat.copy = function copy(A) {
|
592 | return A.map(function(row) {
|
593 | if (isNumber(row))
|
594 | return row;
|
595 | return row.map(function(t) {
|
596 | return t;
|
597 | });
|
598 | });
|
599 | };
|
600 |
|
601 |
|
602 |
|
603 |
|
604 |
|
605 |
|
606 |
|
607 | var jProto = jStat.prototype;
|
608 |
|
609 |
|
610 | jProto.length = 0;
|
611 |
|
612 |
|
613 |
|
614 |
|
615 | jProto.push = Array.prototype.push;
|
616 | jProto.sort = Array.prototype.sort;
|
617 | jProto.splice = Array.prototype.splice;
|
618 | jProto.slice = Array.prototype.slice;
|
619 |
|
620 |
|
621 |
|
622 | jProto.toArray = function toArray() {
|
623 | return this.length > 1 ? slice.call(this) : slice.call(this)[0];
|
624 | };
|
625 |
|
626 |
|
627 |
|
628 | jProto.map = function map(func, toAlter) {
|
629 | return jStat(jStat.map(this, func, toAlter));
|
630 | };
|
631 |
|
632 |
|
633 |
|
634 | jProto.cumreduce = function cumreduce(func, toAlter) {
|
635 | return jStat(jStat.cumreduce(this, func, toAlter));
|
636 | };
|
637 |
|
638 |
|
639 |
|
640 | jProto.alter = function alter(func) {
|
641 | jStat.alter(this, func);
|
642 | return this;
|
643 | };
|
644 |
|
645 |
|
646 |
|
647 | (function(funcs) {
|
648 | for (var i = 0; i < funcs.length; i++) (function(passfunc) {
|
649 | jProto[passfunc] = function(func) {
|
650 | var self = this,
|
651 | results;
|
652 |
|
653 | if (func) {
|
654 | setTimeout(function() {
|
655 | func.call(self, jProto[passfunc].call(self));
|
656 | });
|
657 | return this;
|
658 | }
|
659 | results = jStat[passfunc](this);
|
660 | return isArray(results) ? jStat(results) : results;
|
661 | };
|
662 | })(funcs[i]);
|
663 | })('transpose clear symmetric rows cols dimensions diag antidiag'.split(' '));
|
664 |
|
665 |
|
666 |
|
667 | (function(funcs) {
|
668 | for (var i = 0; i < funcs.length; i++) (function(passfunc) {
|
669 | jProto[passfunc] = function(index, func) {
|
670 | var self = this;
|
671 |
|
672 | if (func) {
|
673 | setTimeout(function() {
|
674 | func.call(self, jProto[passfunc].call(self, index));
|
675 | });
|
676 | return this;
|
677 | }
|
678 | return jStat(jStat[passfunc](this, index));
|
679 | };
|
680 | })(funcs[i]);
|
681 | })('row col'.split(' '));
|
682 |
|
683 |
|
684 |
|
685 | (function(funcs) {
|
686 | for (var i = 0; i < funcs.length; i++) (function(passfunc) {
|
687 | jProto[passfunc] = function() {
|
688 | return jStat(jStat[passfunc].apply(null, arguments));
|
689 | };
|
690 | })(funcs[i]);
|
691 | })('create zeros ones rand identity'.split(' '));
|
692 |
|
693 |
|
694 |
|
695 | return jStat;
|
696 |
|
697 | }(Math));
|
698 | (function(jStat, Math) {
|
699 |
|
700 | var isFunction = jStat.utils.isFunction;
|
701 |
|
702 |
|
703 | function ascNum(a, b) { return a - b; }
|
704 |
|
705 | function clip(arg, min, max) {
|
706 | return Math.max(min, Math.min(arg, max));
|
707 | }
|
708 |
|
709 |
|
710 |
|
711 | jStat.sum = function sum(arr) {
|
712 | var sum = 0;
|
713 | var i = arr.length;
|
714 | while (--i >= 0)
|
715 | sum += arr[i];
|
716 | return sum;
|
717 | };
|
718 |
|
719 |
|
720 |
|
721 | jStat.sumsqrd = function sumsqrd(arr) {
|
722 | var sum = 0;
|
723 | var i = arr.length;
|
724 | while (--i >= 0)
|
725 | sum += arr[i] * arr[i];
|
726 | return sum;
|
727 | };
|
728 |
|
729 |
|
730 |
|
731 | jStat.sumsqerr = function sumsqerr(arr) {
|
732 | var mean = jStat.mean(arr);
|
733 | var sum = 0;
|
734 | var i = arr.length;
|
735 | var tmp;
|
736 | while (--i >= 0) {
|
737 | tmp = arr[i] - mean;
|
738 | sum += tmp * tmp;
|
739 | }
|
740 | return sum;
|
741 | };
|
742 |
|
743 |
|
744 | jStat.sumrow = function sumrow(arr) {
|
745 | var sum = 0;
|
746 | var i = arr.length;
|
747 | while (--i >= 0)
|
748 | sum += arr[i];
|
749 | return sum;
|
750 | };
|
751 |
|
752 |
|
753 | jStat.product = function product(arr) {
|
754 | var prod = 1;
|
755 | var i = arr.length;
|
756 | while (--i >= 0)
|
757 | prod *= arr[i];
|
758 | return prod;
|
759 | };
|
760 |
|
761 |
|
762 |
|
763 | jStat.min = function min(arr) {
|
764 | var low = arr[0];
|
765 | var i = 0;
|
766 | while (++i < arr.length)
|
767 | if (arr[i] < low)
|
768 | low = arr[i];
|
769 | return low;
|
770 | };
|
771 |
|
772 |
|
773 |
|
774 | jStat.max = function max(arr) {
|
775 | var high = arr[0];
|
776 | var i = 0;
|
777 | while (++i < arr.length)
|
778 | if (arr[i] > high)
|
779 | high = arr[i];
|
780 | return high;
|
781 | };
|
782 |
|
783 |
|
784 |
|
785 | jStat.unique = function unique(arr) {
|
786 | var hash = {}, _arr = [];
|
787 | for(var i = 0; i < arr.length; i++) {
|
788 | if (!hash[arr[i]]) {
|
789 | hash[arr[i]] = true;
|
790 | _arr.push(arr[i]);
|
791 | }
|
792 | }
|
793 | return _arr;
|
794 | };
|
795 |
|
796 |
|
797 |
|
798 | jStat.mean = function mean(arr) {
|
799 | return jStat.sum(arr) / arr.length;
|
800 | };
|
801 |
|
802 |
|
803 |
|
804 | jStat.meansqerr = function meansqerr(arr) {
|
805 | return jStat.sumsqerr(arr) / arr.length;
|
806 | };
|
807 |
|
808 |
|
809 |
|
810 | jStat.geomean = function geomean(arr) {
|
811 | return Math.pow(jStat.product(arr), 1 / arr.length);
|
812 | };
|
813 |
|
814 |
|
815 |
|
816 | jStat.median = function median(arr) {
|
817 | var arrlen = arr.length;
|
818 | var _arr = arr.slice().sort(ascNum);
|
819 |
|
820 | return !(arrlen & 1)
|
821 | ? (_arr[(arrlen / 2) - 1 ] + _arr[(arrlen / 2)]) / 2
|
822 | : _arr[(arrlen / 2) | 0 ];
|
823 | };
|
824 |
|
825 |
|
826 |
|
827 | jStat.cumsum = function cumsum(arr) {
|
828 | return jStat.cumreduce(arr, function (a, b) { return a + b; });
|
829 | };
|
830 |
|
831 |
|
832 |
|
833 | jStat.cumprod = function cumprod(arr) {
|
834 | return jStat.cumreduce(arr, function (a, b) { return a * b; });
|
835 | };
|
836 |
|
837 |
|
838 |
|
839 | jStat.diff = function diff(arr) {
|
840 | var diffs = [];
|
841 | var arrLen = arr.length;
|
842 | var i;
|
843 | for (i = 1; i < arrLen; i++)
|
844 | diffs.push(arr[i] - arr[i - 1]);
|
845 | return diffs;
|
846 | };
|
847 |
|
848 |
|
849 |
|
850 | jStat.rank = function (arr) {
|
851 | var arrlen = arr.length;
|
852 | var sorted = arr.slice().sort(ascNum);
|
853 | var ranks = new Array(arrlen);
|
854 | var val;
|
855 | for (var i = 0; i < arrlen; i++) {
|
856 | var first = sorted.indexOf(arr[i]);
|
857 | var last = sorted.lastIndexOf(arr[i]);
|
858 | if (first === last) {
|
859 | val = first;
|
860 | } else {
|
861 | val = (first + last) / 2;
|
862 | }
|
863 | ranks[i] = val + 1;
|
864 | }
|
865 | return ranks;
|
866 | };
|
867 |
|
868 |
|
869 |
|
870 |
|
871 |
|
872 | jStat.mode = function mode(arr) {
|
873 | var arrLen = arr.length;
|
874 | var _arr = arr.slice().sort(ascNum);
|
875 | var count = 1;
|
876 | var maxCount = 0;
|
877 | var numMaxCount = 0;
|
878 | var mode_arr = [];
|
879 | var i;
|
880 |
|
881 | for (i = 0; i < arrLen; i++) {
|
882 | if (_arr[i] === _arr[i + 1]) {
|
883 | count++;
|
884 | } else {
|
885 | if (count > maxCount) {
|
886 | mode_arr = [_arr[i]];
|
887 | maxCount = count;
|
888 | numMaxCount = 0;
|
889 | }
|
890 |
|
891 | else if (count === maxCount) {
|
892 | mode_arr.push(_arr[i]);
|
893 | numMaxCount++;
|
894 | }
|
895 |
|
896 | count = 1;
|
897 | }
|
898 | }
|
899 |
|
900 | return numMaxCount === 0 ? mode_arr[0] : mode_arr;
|
901 | };
|
902 |
|
903 |
|
904 |
|
905 | jStat.range = function range(arr) {
|
906 | return jStat.max(arr) - jStat.min(arr);
|
907 | };
|
908 |
|
909 |
|
910 |
|
911 | jStat.variance = function variance(arr, flag) {
|
912 | return jStat.sumsqerr(arr) / (arr.length - (flag ? 1 : 0));
|
913 | };
|
914 |
|
915 |
|
916 | jStat.pooledvariance = function pooledvariance(arr) {
|
917 | var sumsqerr = arr.reduce(function (a, samples) {return a + jStat.sumsqerr(samples);}, 0);
|
918 | var count = arr.reduce(function (a, samples) {return a + samples.length;}, 0);
|
919 | return sumsqerr / (count - arr.length);
|
920 | };
|
921 |
|
922 |
|
923 | jStat.deviation = function (arr) {
|
924 | var mean = jStat.mean(arr);
|
925 | var arrlen = arr.length;
|
926 | var dev = new Array(arrlen);
|
927 | for (var i = 0; i < arrlen; i++) {
|
928 | dev[i] = arr[i] - mean;
|
929 | }
|
930 | return dev;
|
931 | };
|
932 |
|
933 |
|
934 |
|
935 | jStat.stdev = function stdev(arr, flag) {
|
936 | return Math.sqrt(jStat.variance(arr, flag));
|
937 | };
|
938 |
|
939 |
|
940 | jStat.pooledstdev = function pooledstdev(arr) {
|
941 | return Math.sqrt(jStat.pooledvariance(arr));
|
942 | };
|
943 |
|
944 |
|
945 | jStat.meandev = function meandev(arr) {
|
946 | var mean = jStat.mean(arr);
|
947 | var a = [];
|
948 | for (var i = arr.length - 1; i >= 0; i--) {
|
949 | a.push(Math.abs(arr[i] - mean));
|
950 | }
|
951 | return jStat.mean(a);
|
952 | };
|
953 |
|
954 |
|
955 |
|
956 | jStat.meddev = function meddev(arr) {
|
957 | var median = jStat.median(arr);
|
958 | var a = [];
|
959 | for (var i = arr.length - 1; i >= 0; i--) {
|
960 | a.push(Math.abs(arr[i] - median));
|
961 | }
|
962 | return jStat.median(a);
|
963 | };
|
964 |
|
965 |
|
966 |
|
967 | jStat.coeffvar = function coeffvar(arr) {
|
968 | return jStat.stdev(arr) / jStat.mean(arr);
|
969 | };
|
970 |
|
971 |
|
972 |
|
973 | jStat.quartiles = function quartiles(arr) {
|
974 | var arrlen = arr.length;
|
975 | var _arr = arr.slice().sort(ascNum);
|
976 | return [
|
977 | _arr[ Math.round((arrlen) / 4) - 1 ],
|
978 | _arr[ Math.round((arrlen) / 2) - 1 ],
|
979 | _arr[ Math.round((arrlen) * 3 / 4) - 1 ]
|
980 | ];
|
981 | };
|
982 |
|
983 |
|
984 |
|
985 |
|
986 | jStat.quantiles = function quantiles(arr, quantilesArray, alphap, betap) {
|
987 | var sortedArray = arr.slice().sort(ascNum);
|
988 | var quantileVals = [quantilesArray.length];
|
989 | var n = arr.length;
|
990 | var i, p, m, aleph, k, gamma;
|
991 |
|
992 | if (typeof alphap === 'undefined')
|
993 | alphap = 3 / 8;
|
994 | if (typeof betap === 'undefined')
|
995 | betap = 3 / 8;
|
996 |
|
997 | for (i = 0; i < quantilesArray.length; i++) {
|
998 | p = quantilesArray[i];
|
999 | m = alphap + p * (1 - alphap - betap);
|
1000 | aleph = n * p + m;
|
1001 | k = Math.floor(clip(aleph, 1, n - 1));
|
1002 | gamma = clip(aleph - k, 0, 1);
|
1003 | quantileVals[i] = (1 - gamma) * sortedArray[k - 1] + gamma * sortedArray[k];
|
1004 | }
|
1005 |
|
1006 | return quantileVals;
|
1007 | };
|
1008 |
|
1009 |
|
1010 |
|
1011 | jStat.percentile = function percentile(arr, k, exclusive) {
|
1012 | var _arr = arr.slice().sort(ascNum);
|
1013 | var realIndex = k * (_arr.length + (exclusive ? 1 : -1)) + (exclusive ? 0 : 1);
|
1014 | var index = parseInt(realIndex);
|
1015 | var frac = realIndex - index;
|
1016 | if (index + 1 < _arr.length) {
|
1017 | return _arr[index - 1] + frac * (_arr[index] - _arr[index - 1]);
|
1018 | } else {
|
1019 | return _arr[index - 1];
|
1020 | }
|
1021 | }
|
1022 |
|
1023 |
|
1024 |
|
1025 |
|
1026 | jStat.percentileOfScore = function percentileOfScore(arr, score, kind) {
|
1027 | var counter = 0;
|
1028 | var len = arr.length;
|
1029 | var strict = false;
|
1030 | var value, i;
|
1031 |
|
1032 | if (kind === 'strict')
|
1033 | strict = true;
|
1034 |
|
1035 | for (i = 0; i < len; i++) {
|
1036 | value = arr[i];
|
1037 | if ((strict && value < score) ||
|
1038 | (!strict && value <= score)) {
|
1039 | counter++;
|
1040 | }
|
1041 | }
|
1042 |
|
1043 | return counter / len;
|
1044 | };
|
1045 |
|
1046 |
|
1047 |
|
1048 | jStat.histogram = function histogram(arr, binCnt) {
|
1049 | binCnt = binCnt || 4;
|
1050 | var first = jStat.min(arr);
|
1051 | var binWidth = (jStat.max(arr) - first) / binCnt;
|
1052 | var len = arr.length;
|
1053 | var bins = [];
|
1054 | var i;
|
1055 |
|
1056 | for (i = 0; i < binCnt; i++)
|
1057 | bins[i] = 0;
|
1058 | for (i = 0; i < len; i++)
|
1059 | bins[Math.min(Math.floor(((arr[i] - first) / binWidth)), binCnt - 1)] += 1;
|
1060 |
|
1061 | return bins;
|
1062 | };
|
1063 |
|
1064 |
|
1065 |
|
1066 | jStat.covariance = function covariance(arr1, arr2) {
|
1067 | var u = jStat.mean(arr1);
|
1068 | var v = jStat.mean(arr2);
|
1069 | var arr1Len = arr1.length;
|
1070 | var sq_dev = new Array(arr1Len);
|
1071 | var i;
|
1072 |
|
1073 | for (i = 0; i < arr1Len; i++)
|
1074 | sq_dev[i] = (arr1[i] - u) * (arr2[i] - v);
|
1075 |
|
1076 | return jStat.sum(sq_dev) / (arr1Len - 1);
|
1077 | };
|
1078 |
|
1079 |
|
1080 |
|
1081 | jStat.corrcoeff = function corrcoeff(arr1, arr2) {
|
1082 | return jStat.covariance(arr1, arr2) /
|
1083 | jStat.stdev(arr1, 1) /
|
1084 | jStat.stdev(arr2, 1);
|
1085 | };
|
1086 |
|
1087 |
|
1088 | jStat.spearmancoeff = function (arr1, arr2) {
|
1089 | arr1 = jStat.rank(arr1);
|
1090 | arr2 = jStat.rank(arr2);
|
1091 |
|
1092 | return jStat.corrcoeff(arr1, arr2);
|
1093 | }
|
1094 |
|
1095 |
|
1096 |
|
1097 | jStat.stanMoment = function stanMoment(arr, n) {
|
1098 | var mu = jStat.mean(arr);
|
1099 | var sigma = jStat.stdev(arr);
|
1100 | var len = arr.length;
|
1101 | var skewSum = 0;
|
1102 |
|
1103 | for (var i = 0; i < len; i++)
|
1104 | skewSum += Math.pow((arr[i] - mu) / sigma, n);
|
1105 |
|
1106 | return skewSum / arr.length;
|
1107 | };
|
1108 |
|
1109 |
|
1110 | jStat.skewness = function skewness(arr) {
|
1111 | return jStat.stanMoment(arr, 3);
|
1112 | };
|
1113 |
|
1114 |
|
1115 | jStat.kurtosis = function kurtosis(arr) {
|
1116 | return jStat.stanMoment(arr, 4) - 3;
|
1117 | };
|
1118 |
|
1119 |
|
1120 | var jProto = jStat.prototype;
|
1121 |
|
1122 |
|
1123 |
|
1124 |
|
1125 |
|
1126 |
|
1127 |
|
1128 | (function(funcs) {
|
1129 | for (var i = 0; i < funcs.length; i++) (function(passfunc) {
|
1130 |
|
1131 |
|
1132 | jProto[passfunc] = function(fullbool, func) {
|
1133 | var arr = [];
|
1134 | var i = 0;
|
1135 | var tmpthis = this;
|
1136 |
|
1137 | if (isFunction(fullbool)) {
|
1138 | func = fullbool;
|
1139 | fullbool = false;
|
1140 | }
|
1141 |
|
1142 | if (func) {
|
1143 | setTimeout(function() {
|
1144 | func.call(tmpthis, jProto[passfunc].call(tmpthis, fullbool));
|
1145 | });
|
1146 | return this;
|
1147 | }
|
1148 |
|
1149 | if (this.length > 1) {
|
1150 | tmpthis = fullbool === true ? this : this.transpose();
|
1151 | for (; i < tmpthis.length; i++)
|
1152 | arr[i] = jStat[passfunc](tmpthis[i]);
|
1153 | return arr;
|
1154 | }
|
1155 |
|
1156 | return jStat[passfunc](this[0], fullbool);
|
1157 | };
|
1158 | })(funcs[i]);
|
1159 | })(('cumsum cumprod').split(' '));
|
1160 |
|
1161 |
|
1162 |
|
1163 | (function(funcs) {
|
1164 | for (var i = 0; i < funcs.length; i++) (function(passfunc) {
|
1165 |
|
1166 |
|
1167 | jProto[passfunc] = function(fullbool, func) {
|
1168 | var arr = [];
|
1169 | var i = 0;
|
1170 | var tmpthis = this;
|
1171 |
|
1172 | if (isFunction(fullbool)) {
|
1173 | func = fullbool;
|
1174 | fullbool = false;
|
1175 | }
|
1176 |
|
1177 | if (func) {
|
1178 | setTimeout(function() {
|
1179 | func.call(tmpthis, jProto[passfunc].call(tmpthis, fullbool));
|
1180 | });
|
1181 | return this;
|
1182 | }
|
1183 |
|
1184 | if (this.length > 1) {
|
1185 | if (passfunc !== 'sumrow')
|
1186 | tmpthis = fullbool === true ? this : this.transpose();
|
1187 | for (; i < tmpthis.length; i++)
|
1188 | arr[i] = jStat[passfunc](tmpthis[i]);
|
1189 | return fullbool === true
|
1190 | ? jStat[passfunc](jStat.utils.toVector(arr))
|
1191 | : arr;
|
1192 | }
|
1193 |
|
1194 | return jStat[passfunc](this[0], fullbool);
|
1195 | };
|
1196 | })(funcs[i]);
|
1197 | })(('sum sumsqrd sumsqerr sumrow product min max unique mean meansqerr ' +
|
1198 | 'geomean median diff rank mode range variance deviation stdev meandev ' +
|
1199 | 'meddev coeffvar quartiles histogram skewness kurtosis').split(' '));
|
1200 |
|
1201 |
|
1202 |
|
1203 |
|
1204 | (function(funcs) {
|
1205 | for (var i = 0; i < funcs.length; i++) (function(passfunc) {
|
1206 | jProto[passfunc] = function() {
|
1207 | var arr = [];
|
1208 | var i = 0;
|
1209 | var tmpthis = this;
|
1210 | var args = Array.prototype.slice.call(arguments);
|
1211 | var callbackFunction;
|
1212 |
|
1213 |
|
1214 |
|
1215 | if (isFunction(args[args.length - 1])) {
|
1216 | callbackFunction = args[args.length - 1];
|
1217 | var argsToPass = args.slice(0, args.length - 1);
|
1218 |
|
1219 | setTimeout(function() {
|
1220 | callbackFunction.call(tmpthis,
|
1221 | jProto[passfunc].apply(tmpthis, argsToPass));
|
1222 | });
|
1223 | return this;
|
1224 |
|
1225 |
|
1226 | } else {
|
1227 | callbackFunction = undefined;
|
1228 | var curriedFunction = function curriedFunction(vector) {
|
1229 | return jStat[passfunc].apply(tmpthis, [vector].concat(args));
|
1230 | }
|
1231 | }
|
1232 |
|
1233 |
|
1234 | if (this.length > 1) {
|
1235 | tmpthis = tmpthis.transpose();
|
1236 | for (; i < tmpthis.length; i++)
|
1237 | arr[i] = curriedFunction(tmpthis[i]);
|
1238 | return arr;
|
1239 | }
|
1240 |
|
1241 |
|
1242 | return curriedFunction(this[0]);
|
1243 | };
|
1244 | })(funcs[i]);
|
1245 | })('quantiles percentileOfScore'.split(' '));
|
1246 |
|
1247 | }(jStat, Math));
|
1248 |
|
1249 | (function(jStat, Math) {
|
1250 |
|
1251 |
|
1252 | jStat.gammaln = function gammaln(x) {
|
1253 | var j = 0;
|
1254 | var cof = [
|
1255 | 76.18009172947146, -86.50532032941677, 24.01409824083091,
|
1256 | -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5
|
1257 | ];
|
1258 | var ser = 1.000000000190015;
|
1259 | var xx, y, tmp;
|
1260 | tmp = (y = xx = x) + 5.5;
|
1261 | tmp -= (xx + 0.5) * Math.log(tmp);
|
1262 | for (; j < 6; j++)
|
1263 | ser += cof[j] / ++y;
|
1264 | return Math.log(2.5066282746310005 * ser / xx) - tmp;
|
1265 | };
|
1266 |
|
1267 |
|
1268 |
|
1269 | jStat.gammafn = function gammafn(x) {
|
1270 | var p = [-1.716185138865495, 24.76565080557592, -379.80425647094563,
|
1271 | 629.3311553128184, 866.9662027904133, -31451.272968848367,
|
1272 | -36144.413418691176, 66456.14382024054
|
1273 | ];
|
1274 | var q = [-30.8402300119739, 315.35062697960416, -1015.1563674902192,
|
1275 | -3107.771671572311, 22538.118420980151, 4755.8462775278811,
|
1276 | -134659.9598649693, -115132.2596755535];
|
1277 | var fact = false;
|
1278 | var n = 0;
|
1279 | var xden = 0;
|
1280 | var xnum = 0;
|
1281 | var y = x;
|
1282 | var i, z, yi, res;
|
1283 | if (y <= 0) {
|
1284 | res = y % 1 + 3.6e-16;
|
1285 | if (res) {
|
1286 | fact = (!(y & 1) ? 1 : -1) * Math.PI / Math.sin(Math.PI * res);
|
1287 | y = 1 - y;
|
1288 | } else {
|
1289 | return Infinity;
|
1290 | }
|
1291 | }
|
1292 | yi = y;
|
1293 | if (y < 1) {
|
1294 | z = y++;
|
1295 | } else {
|
1296 | z = (y -= n = (y | 0) - 1) - 1;
|
1297 | }
|
1298 | for (i = 0; i < 8; ++i) {
|
1299 | xnum = (xnum + p[i]) * z;
|
1300 | xden = xden * z + q[i];
|
1301 | }
|
1302 | res = xnum / xden + 1;
|
1303 | if (yi < y) {
|
1304 | res /= yi;
|
1305 | } else if (yi > y) {
|
1306 | for (i = 0; i < n; ++i) {
|
1307 | res *= y;
|
1308 | y++;
|
1309 | }
|
1310 | }
|
1311 | if (fact) {
|
1312 | res = fact / res;
|
1313 | }
|
1314 | return res;
|
1315 | };
|
1316 |
|
1317 |
|
1318 |
|
1319 |
|
1320 | jStat.gammap = function gammap(a, x) {
|
1321 | return jStat.lowRegGamma(a, x) * jStat.gammafn(a);
|
1322 | };
|
1323 |
|
1324 |
|
1325 |
|
1326 | jStat.lowRegGamma = function lowRegGamma(a, x) {
|
1327 | var aln = jStat.gammaln(a);
|
1328 | var ap = a;
|
1329 | var sum = 1 / a;
|
1330 | var del = sum;
|
1331 | var b = x + 1 - a;
|
1332 | var c = 1 / 1.0e-30;
|
1333 | var d = 1 / b;
|
1334 | var h = d;
|
1335 | var i = 1;
|
1336 |
|
1337 | var ITMAX = -~(Math.log((a >= 1) ? a : 1 / a) * 8.5 + a * 0.4 + 17);
|
1338 | var an;
|
1339 |
|
1340 | if (x < 0 || a <= 0) {
|
1341 | return NaN;
|
1342 | } else if (x < a + 1) {
|
1343 | for (; i <= ITMAX; i++) {
|
1344 | sum += del *= x / ++ap;
|
1345 | }
|
1346 | return (sum * Math.exp(-x + a * Math.log(x) - (aln)));
|
1347 | }
|
1348 |
|
1349 | for (; i <= ITMAX; i++) {
|
1350 | an = -i * (i - a);
|
1351 | b += 2;
|
1352 | d = an * d + b;
|
1353 | c = b + an / c;
|
1354 | d = 1 / d;
|
1355 | h *= d * c;
|
1356 | }
|
1357 |
|
1358 | return (1 - h * Math.exp(-x + a * Math.log(x) - (aln)));
|
1359 | };
|
1360 |
|
1361 |
|
1362 | jStat.factorialln = function factorialln(n) {
|
1363 | return n < 0 ? NaN : jStat.gammaln(n + 1);
|
1364 | };
|
1365 |
|
1366 |
|
1367 | jStat.factorial = function factorial(n) {
|
1368 | return n < 0 ? NaN : jStat.gammafn(n + 1);
|
1369 | };
|
1370 |
|
1371 |
|
1372 | jStat.combination = function combination(n, m) {
|
1373 |
|
1374 | return (n > 170 || m > 170)
|
1375 | ? Math.exp(jStat.combinationln(n, m))
|
1376 | : (jStat.factorial(n) / jStat.factorial(m)) / jStat.factorial(n - m);
|
1377 | };
|
1378 |
|
1379 |
|
1380 | jStat.combinationln = function combinationln(n, m){
|
1381 | return jStat.factorialln(n) - jStat.factorialln(m) - jStat.factorialln(n - m);
|
1382 | };
|
1383 |
|
1384 |
|
1385 |
|
1386 | jStat.permutation = function permutation(n, m) {
|
1387 | return jStat.factorial(n) / jStat.factorial(n - m);
|
1388 | };
|
1389 |
|
1390 |
|
1391 |
|
1392 | jStat.betafn = function betafn(x, y) {
|
1393 |
|
1394 | if (x <= 0 || y <= 0)
|
1395 | return undefined;
|
1396 |
|
1397 | return (x + y > 170)
|
1398 | ? Math.exp(jStat.betaln(x, y))
|
1399 | : jStat.gammafn(x) * jStat.gammafn(y) / jStat.gammafn(x + y);
|
1400 | };
|
1401 |
|
1402 |
|
1403 |
|
1404 | jStat.betaln = function betaln(x, y) {
|
1405 | return jStat.gammaln(x) + jStat.gammaln(y) - jStat.gammaln(x + y);
|
1406 | };
|
1407 |
|
1408 |
|
1409 |
|
1410 |
|
1411 | jStat.betacf = function betacf(x, a, b) {
|
1412 | var fpmin = 1e-30;
|
1413 | var m = 1;
|
1414 | var qab = a + b;
|
1415 | var qap = a + 1;
|
1416 | var qam = a - 1;
|
1417 | var c = 1;
|
1418 | var d = 1 - qab * x / qap;
|
1419 | var m2, aa, del, h;
|
1420 |
|
1421 |
|
1422 | if (Math.abs(d) < fpmin)
|
1423 | d = fpmin;
|
1424 | d = 1 / d;
|
1425 | h = d;
|
1426 |
|
1427 | for (; m <= 100; m++) {
|
1428 | m2 = 2 * m;
|
1429 | aa = m * (b - m) * x / ((qam + m2) * (a + m2));
|
1430 |
|
1431 | d = 1 + aa * d;
|
1432 | if (Math.abs(d) < fpmin)
|
1433 | d = fpmin;
|
1434 | c = 1 + aa / c;
|
1435 | if (Math.abs(c) < fpmin)
|
1436 | c = fpmin;
|
1437 | d = 1 / d;
|
1438 | h *= d * c;
|
1439 | aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
|
1440 |
|
1441 | d = 1 + aa * d;
|
1442 | if (Math.abs(d) < fpmin)
|
1443 | d = fpmin;
|
1444 | c = 1 + aa / c;
|
1445 | if (Math.abs(c) < fpmin)
|
1446 | c = fpmin;
|
1447 | d = 1 / d;
|
1448 | del = d * c;
|
1449 | h *= del;
|
1450 | if (Math.abs(del - 1.0) < 3e-7)
|
1451 | break;
|
1452 | }
|
1453 |
|
1454 | return h;
|
1455 | };
|
1456 |
|
1457 |
|
1458 |
|
1459 | jStat.gammapinv = function gammapinv(p, a) {
|
1460 | var j = 0;
|
1461 | var a1 = a - 1;
|
1462 | var EPS = 1e-8;
|
1463 | var gln = jStat.gammaln(a);
|
1464 | var x, err, t, u, pp, lna1, afac;
|
1465 |
|
1466 | if (p >= 1)
|
1467 | return Math.max(100, a + 100 * Math.sqrt(a));
|
1468 | if (p <= 0)
|
1469 | return 0;
|
1470 | if (a > 1) {
|
1471 | lna1 = Math.log(a1);
|
1472 | afac = Math.exp(a1 * (lna1 - 1) - gln);
|
1473 | pp = (p < 0.5) ? p : 1 - p;
|
1474 | t = Math.sqrt(-2 * Math.log(pp));
|
1475 | x = (2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t;
|
1476 | if (p < 0.5)
|
1477 | x = -x;
|
1478 | x = Math.max(1e-3,
|
1479 | a * Math.pow(1 - 1 / (9 * a) - x / (3 * Math.sqrt(a)), 3));
|
1480 | } else {
|
1481 | t = 1 - a * (0.253 + a * 0.12);
|
1482 | if (p < t)
|
1483 | x = Math.pow(p / t, 1 / a);
|
1484 | else
|
1485 | x = 1 - Math.log(1 - (p - t) / (1 - t));
|
1486 | }
|
1487 |
|
1488 | for(; j < 12; j++) {
|
1489 | if (x <= 0)
|
1490 | return 0;
|
1491 | err = jStat.lowRegGamma(a, x) - p;
|
1492 | if (a > 1)
|
1493 | t = afac * Math.exp(-(x - a1) + a1 * (Math.log(x) - lna1));
|
1494 | else
|
1495 | t = Math.exp(-x + a1 * Math.log(x) - gln);
|
1496 | u = err / t;
|
1497 | x -= (t = u / (1 - 0.5 * Math.min(1, u * ((a - 1) / x - 1))));
|
1498 | if (x <= 0)
|
1499 | x = 0.5 * (x + t);
|
1500 | if (Math.abs(t) < EPS * x)
|
1501 | break;
|
1502 | }
|
1503 |
|
1504 | return x;
|
1505 | };
|
1506 |
|
1507 |
|
1508 |
|
1509 | jStat.erf = function erf(x) {
|
1510 | var cof = [-1.3026537197817094, 6.4196979235649026e-1, 1.9476473204185836e-2,
|
1511 | -9.561514786808631e-3, -9.46595344482036e-4, 3.66839497852761e-4,
|
1512 | 4.2523324806907e-5, -2.0278578112534e-5, -1.624290004647e-6,
|
1513 | 1.303655835580e-6, 1.5626441722e-8, -8.5238095915e-8,
|
1514 | 6.529054439e-9, 5.059343495e-9, -9.91364156e-10,
|
1515 | -2.27365122e-10, 9.6467911e-11, 2.394038e-12,
|
1516 | -6.886027e-12, 8.94487e-13, 3.13092e-13,
|
1517 | -1.12708e-13, 3.81e-16, 7.106e-15,
|
1518 | -1.523e-15, -9.4e-17, 1.21e-16,
|
1519 | -2.8e-17];
|
1520 | var j = cof.length - 1;
|
1521 | var isneg = false;
|
1522 | var d = 0;
|
1523 | var dd = 0;
|
1524 | var t, ty, tmp, res;
|
1525 |
|
1526 | if (x < 0) {
|
1527 | x = -x;
|
1528 | isneg = true;
|
1529 | }
|
1530 |
|
1531 | t = 2 / (2 + x);
|
1532 | ty = 4 * t - 2;
|
1533 |
|
1534 | for(; j > 0; j--) {
|
1535 | tmp = d;
|
1536 | d = ty * d - dd + cof[j];
|
1537 | dd = tmp;
|
1538 | }
|
1539 |
|
1540 | res = t * Math.exp(-x * x + 0.5 * (cof[0] + ty * d) - dd);
|
1541 | return isneg ? res - 1 : 1 - res;
|
1542 | };
|
1543 |
|
1544 |
|
1545 |
|
1546 | jStat.erfc = function erfc(x) {
|
1547 | return 1 - jStat.erf(x);
|
1548 | };
|
1549 |
|
1550 |
|
1551 |
|
1552 | jStat.erfcinv = function erfcinv(p) {
|
1553 | var j = 0;
|
1554 | var x, err, t, pp;
|
1555 | if (p >= 2)
|
1556 | return -100;
|
1557 | if (p <= 0)
|
1558 | return 100;
|
1559 | pp = (p < 1) ? p : 2 - p;
|
1560 | t = Math.sqrt(-2 * Math.log(pp / 2));
|
1561 | x = -0.70711 * ((2.30753 + t * 0.27061) /
|
1562 | (1 + t * (0.99229 + t * 0.04481)) - t);
|
1563 | for (; j < 2; j++) {
|
1564 | err = jStat.erfc(x) - pp;
|
1565 | x += err / (1.12837916709551257 * Math.exp(-x * x) - x * err);
|
1566 | }
|
1567 | return (p < 1) ? x : -x;
|
1568 | };
|
1569 |
|
1570 |
|
1571 |
|
1572 | jStat.ibetainv = function ibetainv(p, a, b) {
|
1573 | var EPS = 1e-8;
|
1574 | var a1 = a - 1;
|
1575 | var b1 = b - 1;
|
1576 | var j = 0;
|
1577 | var lna, lnb, pp, t, u, err, x, al, h, w, afac;
|
1578 | if (p <= 0)
|
1579 | return 0;
|
1580 | if (p >= 1)
|
1581 | return 1;
|
1582 | if (a >= 1 && b >= 1) {
|
1583 | pp = (p < 0.5) ? p : 1 - p;
|
1584 | t = Math.sqrt(-2 * Math.log(pp));
|
1585 | x = (2.30753 + t * 0.27061) / (1 + t* (0.99229 + t * 0.04481)) - t;
|
1586 | if (p < 0.5)
|
1587 | x = -x;
|
1588 | al = (x * x - 3) / 6;
|
1589 | h = 2 / (1 / (2 * a - 1) + 1 / (2 * b - 1));
|
1590 | w = (x * Math.sqrt(al + h) / h) - (1 / (2 * b - 1) - 1 / (2 * a - 1)) *
|
1591 | (al + 5 / 6 - 2 / (3 * h));
|
1592 | x = a / (a + b * Math.exp(2 * w));
|
1593 | } else {
|
1594 | lna = Math.log(a / (a + b));
|
1595 | lnb = Math.log(b / (a + b));
|
1596 | t = Math.exp(a * lna) / a;
|
1597 | u = Math.exp(b * lnb) / b;
|
1598 | w = t + u;
|
1599 | if (p < t / w)
|
1600 | x = Math.pow(a * w * p, 1 / a);
|
1601 | else
|
1602 | x = 1 - Math.pow(b * w * (1 - p), 1 / b);
|
1603 | }
|
1604 | afac = -jStat.gammaln(a) - jStat.gammaln(b) + jStat.gammaln(a + b);
|
1605 | for(; j < 10; j++) {
|
1606 | if (x === 0 || x === 1)
|
1607 | return x;
|
1608 | err = jStat.ibeta(x, a, b) - p;
|
1609 | t = Math.exp(a1 * Math.log(x) + b1 * Math.log(1 - x) + afac);
|
1610 | u = err / t;
|
1611 | x -= (t = u / (1 - 0.5 * Math.min(1, u * (a1 / x - b1 / (1 - x)))));
|
1612 | if (x <= 0)
|
1613 | x = 0.5 * (x + t);
|
1614 | if (x >= 1)
|
1615 | x = 0.5 * (x + t + 1);
|
1616 | if (Math.abs(t) < EPS * x && j > 0)
|
1617 | break;
|
1618 | }
|
1619 | return x;
|
1620 | };
|
1621 |
|
1622 |
|
1623 |
|
1624 | jStat.ibeta = function ibeta(x, a, b) {
|
1625 |
|
1626 | var bt = (x === 0 || x === 1) ? 0 :
|
1627 | Math.exp(jStat.gammaln(a + b) - jStat.gammaln(a) -
|
1628 | jStat.gammaln(b) + a * Math.log(x) + b *
|
1629 | Math.log(1 - x));
|
1630 | if (x < 0 || x > 1)
|
1631 | return false;
|
1632 | if (x < (a + 1) / (a + b + 2))
|
1633 |
|
1634 | return bt * jStat.betacf(x, a, b) / a;
|
1635 |
|
1636 | return 1 - bt * jStat.betacf(1 - x, b, a) / b;
|
1637 | };
|
1638 |
|
1639 |
|
1640 |
|
1641 |
|
1642 | jStat.randn = function randn(n, m) {
|
1643 | var u, v, x, y, q;
|
1644 | if (!m)
|
1645 | m = n;
|
1646 | if (n)
|
1647 | return jStat.create(n, m, function() { return jStat.randn(); });
|
1648 | do {
|
1649 | u = jStat._random_fn();
|
1650 | v = 1.7156 * (jStat._random_fn() - 0.5);
|
1651 | x = u - 0.449871;
|
1652 | y = Math.abs(v) + 0.386595;
|
1653 | q = x * x + y * (0.19600 * y - 0.25472 * x);
|
1654 | } while (q > 0.27597 && (q > 0.27846 || v * v > -4 * Math.log(u) * u * u));
|
1655 | return v / u;
|
1656 | };
|
1657 |
|
1658 |
|
1659 |
|
1660 | jStat.randg = function randg(shape, n, m) {
|
1661 | var oalph = shape;
|
1662 | var a1, a2, u, v, x, mat;
|
1663 | if (!m)
|
1664 | m = n;
|
1665 | if (!shape)
|
1666 | shape = 1;
|
1667 | if (n) {
|
1668 | mat = jStat.zeros(n,m);
|
1669 | mat.alter(function() { return jStat.randg(shape); });
|
1670 | return mat;
|
1671 | }
|
1672 | if (shape < 1)
|
1673 | shape += 1;
|
1674 | a1 = shape - 1 / 3;
|
1675 | a2 = 1 / Math.sqrt(9 * a1);
|
1676 | do {
|
1677 | do {
|
1678 | x = jStat.randn();
|
1679 | v = 1 + a2 * x;
|
1680 | } while(v <= 0);
|
1681 | v = v * v * v;
|
1682 | u = jStat._random_fn();
|
1683 | } while(u > 1 - 0.331 * Math.pow(x, 4) &&
|
1684 | Math.log(u) > 0.5 * x*x + a1 * (1 - v + Math.log(v)));
|
1685 |
|
1686 | if (shape == oalph)
|
1687 | return a1 * v;
|
1688 |
|
1689 | do {
|
1690 | u = jStat._random_fn();
|
1691 | } while(u === 0);
|
1692 | return Math.pow(u, 1 / oalph) * a1 * v;
|
1693 | };
|
1694 |
|
1695 |
|
1696 |
|
1697 | (function(funcs) {
|
1698 | for (var i = 0; i < funcs.length; i++) (function(passfunc) {
|
1699 | jStat.fn[passfunc] = function() {
|
1700 | return jStat(
|
1701 | jStat.map(this, function(value) { return jStat[passfunc](value); }));
|
1702 | }
|
1703 | })(funcs[i]);
|
1704 | })('gammaln gammafn factorial factorialln'.split(' '));
|
1705 |
|
1706 |
|
1707 | (function(funcs) {
|
1708 | for (var i = 0; i < funcs.length; i++) (function(passfunc) {
|
1709 | jStat.fn[passfunc] = function() {
|
1710 | return jStat(jStat[passfunc].apply(null, arguments));
|
1711 | };
|
1712 | })(funcs[i]);
|
1713 | })('randn'.split(' '));
|
1714 |
|
1715 | }(jStat, Math));
|
1716 | (function(jStat, Math) {
|
1717 |
|
1718 |
|
1719 | (function(list) {
|
1720 | for (var i = 0; i < list.length; i++) (function(func) {
|
1721 |
|
1722 | jStat[func] = function(a, b, c) {
|
1723 | if (!(this instanceof arguments.callee))
|
1724 | return new arguments.callee(a, b, c);
|
1725 | this._a = a;
|
1726 | this._b = b;
|
1727 | this._c = c;
|
1728 | return this;
|
1729 | };
|
1730 |
|
1731 | jStat.fn[func] = function(a, b, c) {
|
1732 | var newthis = jStat[func](a, b, c);
|
1733 | newthis.data = this;
|
1734 | return newthis;
|
1735 | };
|
1736 |
|
1737 | jStat[func].prototype.sample = function(arr) {
|
1738 | var a = this._a;
|
1739 | var b = this._b;
|
1740 | var c = this._c;
|
1741 | if (arr)
|
1742 | return jStat.alter(arr, function() {
|
1743 | return jStat[func].sample(a, b, c);
|
1744 | });
|
1745 | else
|
1746 | return jStat[func].sample(a, b, c);
|
1747 | };
|
1748 |
|
1749 | (function(vals) {
|
1750 | for (var i = 0; i < vals.length; i++) (function(fnfunc) {
|
1751 | jStat[func].prototype[fnfunc] = function(x) {
|
1752 | var a = this._a;
|
1753 | var b = this._b;
|
1754 | var c = this._c;
|
1755 | if (!x && x !== 0)
|
1756 | x = this.data;
|
1757 | if (typeof x !== 'number') {
|
1758 | return jStat.fn.map.call(x, function(x) {
|
1759 | return jStat[func][fnfunc](x, a, b, c);
|
1760 | });
|
1761 | }
|
1762 | return jStat[func][fnfunc](x, a, b, c);
|
1763 | };
|
1764 | })(vals[i]);
|
1765 | })('pdf cdf inv'.split(' '));
|
1766 |
|
1767 | (function(vals) {
|
1768 | for (var i = 0; i < vals.length; i++) (function(fnfunc) {
|
1769 | jStat[func].prototype[fnfunc] = function() {
|
1770 | return jStat[func][fnfunc](this._a, this._b, this._c);
|
1771 | };
|
1772 | })(vals[i]);
|
1773 | })('mean median mode variance'.split(' '));
|
1774 | })(list[i]);
|
1775 | })((
|
1776 | 'beta centralF cauchy chisquare exponential gamma invgamma kumaraswamy ' +
|
1777 | 'laplace lognormal noncentralt normal pareto studentt weibull uniform ' +
|
1778 | 'binomial negbin hypgeom poisson triangular tukey arcsine'
|
1779 | ).split(' '));
|
1780 |
|
1781 |
|
1782 |
|
1783 |
|
1784 | jStat.extend(jStat.beta, {
|
1785 | pdf: function pdf(x, alpha, beta) {
|
1786 |
|
1787 | if (x > 1 || x < 0)
|
1788 | return 0;
|
1789 |
|
1790 | if (alpha == 1 && beta == 1)
|
1791 | return 1;
|
1792 |
|
1793 | if (alpha < 512 && beta < 512) {
|
1794 | return (Math.pow(x, alpha - 1) * Math.pow(1 - x, beta - 1)) /
|
1795 | jStat.betafn(alpha, beta);
|
1796 | } else {
|
1797 | return Math.exp((alpha - 1) * Math.log(x) +
|
1798 | (beta - 1) * Math.log(1 - x) -
|
1799 | jStat.betaln(alpha, beta));
|
1800 | }
|
1801 | },
|
1802 |
|
1803 | cdf: function cdf(x, alpha, beta) {
|
1804 | return (x > 1 || x < 0) ? (x > 1) * 1 : jStat.ibeta(x, alpha, beta);
|
1805 | },
|
1806 |
|
1807 | inv: function inv(x, alpha, beta) {
|
1808 | return jStat.ibetainv(x, alpha, beta);
|
1809 | },
|
1810 |
|
1811 | mean: function mean(alpha, beta) {
|
1812 | return alpha / (alpha + beta);
|
1813 | },
|
1814 |
|
1815 | median: function median(alpha, beta) {
|
1816 | return jStat.ibetainv(0.5, alpha, beta);
|
1817 | },
|
1818 |
|
1819 | mode: function mode(alpha, beta) {
|
1820 | return (alpha - 1 ) / ( alpha + beta - 2);
|
1821 | },
|
1822 |
|
1823 |
|
1824 | sample: function sample(alpha, beta) {
|
1825 | var u = jStat.randg(alpha);
|
1826 | return u / (u + jStat.randg(beta));
|
1827 | },
|
1828 |
|
1829 | variance: function variance(alpha, beta) {
|
1830 | return (alpha * beta) / (Math.pow(alpha + beta, 2) * (alpha + beta + 1));
|
1831 | }
|
1832 | });
|
1833 |
|
1834 |
|
1835 | jStat.extend(jStat.centralF, {
|
1836 |
|
1837 |
|
1838 |
|
1839 | pdf: function pdf(x, df1, df2) {
|
1840 | var p, q, f;
|
1841 |
|
1842 | if (x < 0)
|
1843 | return 0;
|
1844 |
|
1845 | if (df1 <= 2) {
|
1846 | if (x === 0 && df1 < 2) {
|
1847 | return Infinity;
|
1848 | }
|
1849 | if (x === 0 && df1 === 2) {
|
1850 | return 1;
|
1851 | }
|
1852 | return (1 / jStat.betafn(df1 / 2, df2 / 2)) *
|
1853 | Math.pow(df1 / df2, df1 / 2) *
|
1854 | Math.pow(x, (df1/2) - 1) *
|
1855 | Math.pow((1 + (df1 / df2) * x), -(df1 + df2) / 2);
|
1856 | }
|
1857 |
|
1858 | p = (df1 * x) / (df2 + x * df1);
|
1859 | q = df2 / (df2 + x * df1);
|
1860 | f = df1 * q / 2.0;
|
1861 | return f * jStat.binomial.pdf((df1 - 2) / 2, (df1 + df2 - 2) / 2, p);
|
1862 | },
|
1863 |
|
1864 | cdf: function cdf(x, df1, df2) {
|
1865 | if (x < 0)
|
1866 | return 0;
|
1867 | return jStat.ibeta((df1 * x) / (df1 * x + df2), df1 / 2, df2 / 2);
|
1868 | },
|
1869 |
|
1870 | inv: function inv(x, df1, df2) {
|
1871 | return df2 / (df1 * (1 / jStat.ibetainv(x, df1 / 2, df2 / 2) - 1));
|
1872 | },
|
1873 |
|
1874 | mean: function mean(df1, df2) {
|
1875 | return (df2 > 2) ? df2 / (df2 - 2) : undefined;
|
1876 | },
|
1877 |
|
1878 | mode: function mode(df1, df2) {
|
1879 | return (df1 > 2) ? (df2 * (df1 - 2)) / (df1 * (df2 + 2)) : undefined;
|
1880 | },
|
1881 |
|
1882 |
|
1883 | sample: function sample(df1, df2) {
|
1884 | var x1 = jStat.randg(df1 / 2) * 2;
|
1885 | var x2 = jStat.randg(df2 / 2) * 2;
|
1886 | return (x1 / df1) / (x2 / df2);
|
1887 | },
|
1888 |
|
1889 | variance: function variance(df1, df2) {
|
1890 | if (df2 <= 4)
|
1891 | return undefined;
|
1892 | return 2 * df2 * df2 * (df1 + df2 - 2) /
|
1893 | (df1 * (df2 - 2) * (df2 - 2) * (df2 - 4));
|
1894 | }
|
1895 | });
|
1896 |
|
1897 |
|
1898 |
|
1899 | jStat.extend(jStat.cauchy, {
|
1900 | pdf: function pdf(x, local, scale) {
|
1901 | if (scale < 0) { return 0; }
|
1902 |
|
1903 | return (scale / (Math.pow(x - local, 2) + Math.pow(scale, 2))) / Math.PI;
|
1904 | },
|
1905 |
|
1906 | cdf: function cdf(x, local, scale) {
|
1907 | return Math.atan((x - local) / scale) / Math.PI + 0.5;
|
1908 | },
|
1909 |
|
1910 | inv: function(p, local, scale) {
|
1911 | return local + scale * Math.tan(Math.PI * (p - 0.5));
|
1912 | },
|
1913 |
|
1914 | median: function median(local/*, scale*/) {
|
1915 | return local;
|
1916 | },
|
1917 |
|
1918 | mode: function mode(local/*, scale*/) {
|
1919 | return local;
|
1920 | },
|
1921 |
|
1922 | sample: function sample(local, scale) {
|
1923 | return jStat.randn() *
|
1924 | Math.sqrt(1 / (2 * jStat.randg(0.5))) * scale + local;
|
1925 | }
|
1926 | });
|
1927 |
|
1928 |
|
1929 |
|
1930 |
|
1931 | jStat.extend(jStat.chisquare, {
|
1932 | pdf: function pdf(x, dof) {
|
1933 | if (x < 0)
|
1934 | return 0;
|
1935 | return (x === 0 && dof === 2) ? 0.5 :
|
1936 | Math.exp((dof / 2 - 1) * Math.log(x) - x / 2 - (dof / 2) *
|
1937 | Math.log(2) - jStat.gammaln(dof / 2));
|
1938 | },
|
1939 |
|
1940 | cdf: function cdf(x, dof) {
|
1941 | if (x < 0)
|
1942 | return 0;
|
1943 | return jStat.lowRegGamma(dof / 2, x / 2);
|
1944 | },
|
1945 |
|
1946 | inv: function(p, dof) {
|
1947 | return 2 * jStat.gammapinv(p, 0.5 * dof);
|
1948 | },
|
1949 |
|
1950 | mean : function(dof) {
|
1951 | return dof;
|
1952 | },
|
1953 |
|
1954 |
|
1955 | median: function median(dof) {
|
1956 | return dof * Math.pow(1 - (2 / (9 * dof)), 3);
|
1957 | },
|
1958 |
|
1959 | mode: function mode(dof) {
|
1960 | return (dof - 2 > 0) ? dof - 2 : 0;
|
1961 | },
|
1962 |
|
1963 | sample: function sample(dof) {
|
1964 | return jStat.randg(dof / 2) * 2;
|
1965 | },
|
1966 |
|
1967 | variance: function variance(dof) {
|
1968 | return 2 * dof;
|
1969 | }
|
1970 | });
|
1971 |
|
1972 |
|
1973 |
|
1974 |
|
1975 | jStat.extend(jStat.exponential, {
|
1976 | pdf: function pdf(x, rate) {
|
1977 | return x < 0 ? 0 : rate * Math.exp(-rate * x);
|
1978 | },
|
1979 |
|
1980 | cdf: function cdf(x, rate) {
|
1981 | return x < 0 ? 0 : 1 - Math.exp(-rate * x);
|
1982 | },
|
1983 |
|
1984 | inv: function(p, rate) {
|
1985 | return -Math.log(1 - p) / rate;
|
1986 | },
|
1987 |
|
1988 | mean : function(rate) {
|
1989 | return 1 / rate;
|
1990 | },
|
1991 |
|
1992 | median: function (rate) {
|
1993 | return (1 / rate) * Math.log(2);
|
1994 | },
|
1995 |
|
1996 | mode: function mode(/*rate*/) {
|
1997 | return 0;
|
1998 | },
|
1999 |
|
2000 | sample: function sample(rate) {
|
2001 | return -1 / rate * Math.log(jStat._random_fn());
|
2002 | },
|
2003 |
|
2004 | variance : function(rate) {
|
2005 | return Math.pow(rate, -2);
|
2006 | }
|
2007 | });
|
2008 |
|
2009 |
|
2010 |
|
2011 |
|
2012 | jStat.extend(jStat.gamma, {
|
2013 | pdf: function pdf(x, shape, scale) {
|
2014 | if (x < 0)
|
2015 | return 0;
|
2016 | return (x === 0 && shape === 1) ? 1 / scale :
|
2017 | Math.exp((shape - 1) * Math.log(x) - x / scale -
|
2018 | jStat.gammaln(shape) - shape * Math.log(scale));
|
2019 | },
|
2020 |
|
2021 | cdf: function cdf(x, shape, scale) {
|
2022 | if (x < 0)
|
2023 | return 0;
|
2024 | return jStat.lowRegGamma(shape, x / scale);
|
2025 | },
|
2026 |
|
2027 | inv: function(p, shape, scale) {
|
2028 | return jStat.gammapinv(p, shape) * scale;
|
2029 | },
|
2030 |
|
2031 | mean : function(shape, scale) {
|
2032 | return shape * scale;
|
2033 | },
|
2034 |
|
2035 | mode: function mode(shape, scale) {
|
2036 | if(shape > 1) return (shape - 1) * scale;
|
2037 | return undefined;
|
2038 | },
|
2039 |
|
2040 | sample: function sample(shape, scale) {
|
2041 | return jStat.randg(shape) * scale;
|
2042 | },
|
2043 |
|
2044 | variance: function variance(shape, scale) {
|
2045 | return shape * scale * scale;
|
2046 | }
|
2047 | });
|
2048 |
|
2049 |
|
2050 | jStat.extend(jStat.invgamma, {
|
2051 | pdf: function pdf(x, shape, scale) {
|
2052 | if (x <= 0)
|
2053 | return 0;
|
2054 | return Math.exp(-(shape + 1) * Math.log(x) - scale / x -
|
2055 | jStat.gammaln(shape) + shape * Math.log(scale));
|
2056 | },
|
2057 |
|
2058 | cdf: function cdf(x, shape, scale) {
|
2059 | if (x <= 0)
|
2060 | return 0;
|
2061 | return 1 - jStat.lowRegGamma(shape, scale / x);
|
2062 | },
|
2063 |
|
2064 | inv: function(p, shape, scale) {
|
2065 | return scale / jStat.gammapinv(1 - p, shape);
|
2066 | },
|
2067 |
|
2068 | mean : function(shape, scale) {
|
2069 | return (shape > 1) ? scale / (shape - 1) : undefined;
|
2070 | },
|
2071 |
|
2072 | mode: function mode(shape, scale) {
|
2073 | return scale / (shape + 1);
|
2074 | },
|
2075 |
|
2076 | sample: function sample(shape, scale) {
|
2077 | return scale / jStat.randg(shape);
|
2078 | },
|
2079 |
|
2080 | variance: function variance(shape, scale) {
|
2081 | if (shape <= 2)
|
2082 | return undefined;
|
2083 | return scale * scale / ((shape - 1) * (shape - 1) * (shape - 2));
|
2084 | }
|
2085 | });
|
2086 |
|
2087 |
|
2088 |
|
2089 | jStat.extend(jStat.kumaraswamy, {
|
2090 | pdf: function pdf(x, alpha, beta) {
|
2091 | if (x === 0 && alpha === 1)
|
2092 | return beta;
|
2093 | else if (x === 1 && beta === 1)
|
2094 | return alpha;
|
2095 | return Math.exp(Math.log(alpha) + Math.log(beta) + (alpha - 1) *
|
2096 | Math.log(x) + (beta - 1) *
|
2097 | Math.log(1 - Math.pow(x, alpha)));
|
2098 | },
|
2099 |
|
2100 | cdf: function cdf(x, alpha, beta) {
|
2101 | if (x < 0)
|
2102 | return 0;
|
2103 | else if (x > 1)
|
2104 | return 1;
|
2105 | return (1 - Math.pow(1 - Math.pow(x, alpha), beta));
|
2106 | },
|
2107 |
|
2108 | inv: function inv(p, alpha, beta) {
|
2109 | return Math.pow(1 - Math.pow(1 - p, 1 / beta), 1 / alpha);
|
2110 | },
|
2111 |
|
2112 | mean : function(alpha, beta) {
|
2113 | return (beta * jStat.gammafn(1 + 1 / alpha) *
|
2114 | jStat.gammafn(beta)) / (jStat.gammafn(1 + 1 / alpha + beta));
|
2115 | },
|
2116 |
|
2117 | median: function median(alpha, beta) {
|
2118 | return Math.pow(1 - Math.pow(2, -1 / beta), 1 / alpha);
|
2119 | },
|
2120 |
|
2121 | mode: function mode(alpha, beta) {
|
2122 | if (!(alpha >= 1 && beta >= 1 && (alpha !== 1 && beta !== 1)))
|
2123 | return undefined;
|
2124 | return Math.pow((alpha - 1) / (alpha * beta - 1), 1 / alpha);
|
2125 | },
|
2126 |
|
2127 | variance: function variance(/*alpha, beta*/) {
|
2128 | throw new Error('variance not yet implemented');
|
2129 |
|
2130 | }
|
2131 | });
|
2132 |
|
2133 |
|
2134 |
|
2135 |
|
2136 | jStat.extend(jStat.lognormal, {
|
2137 | pdf: function pdf(x, mu, sigma) {
|
2138 | if (x <= 0)
|
2139 | return 0;
|
2140 | return Math.exp(-Math.log(x) - 0.5 * Math.log(2 * Math.PI) -
|
2141 | Math.log(sigma) - Math.pow(Math.log(x) - mu, 2) /
|
2142 | (2 * sigma * sigma));
|
2143 | },
|
2144 |
|
2145 | cdf: function cdf(x, mu, sigma) {
|
2146 | if (x < 0)
|
2147 | return 0;
|
2148 | return 0.5 +
|
2149 | (0.5 * jStat.erf((Math.log(x) - mu) / Math.sqrt(2 * sigma * sigma)));
|
2150 | },
|
2151 |
|
2152 | inv: function(p, mu, sigma) {
|
2153 | return Math.exp(-1.41421356237309505 * sigma * jStat.erfcinv(2 * p) + mu);
|
2154 | },
|
2155 |
|
2156 | mean: function mean(mu, sigma) {
|
2157 | return Math.exp(mu + sigma * sigma / 2);
|
2158 | },
|
2159 |
|
2160 | median: function median(mu/*, sigma*/) {
|
2161 | return Math.exp(mu);
|
2162 | },
|
2163 |
|
2164 | mode: function mode(mu, sigma) {
|
2165 | return Math.exp(mu - sigma * sigma);
|
2166 | },
|
2167 |
|
2168 | sample: function sample(mu, sigma) {
|
2169 | return Math.exp(jStat.randn() * sigma + mu);
|
2170 | },
|
2171 |
|
2172 | variance: function variance(mu, sigma) {
|
2173 | return (Math.exp(sigma * sigma) - 1) * Math.exp(2 * mu + sigma * sigma);
|
2174 | }
|
2175 | });
|
2176 |
|
2177 |
|
2178 |
|
2179 |
|
2180 | jStat.extend(jStat.noncentralt, {
|
2181 | pdf: function pdf(x, dof, ncp) {
|
2182 | var tol = 1e-14;
|
2183 | if (Math.abs(ncp) < tol)
|
2184 | return jStat.studentt.pdf(x, dof)
|
2185 |
|
2186 | if (Math.abs(x) < tol) {
|
2187 | return Math.exp(jStat.gammaln((dof + 1) / 2) - ncp * ncp / 2 -
|
2188 | 0.5 * Math.log(Math.PI * dof) - jStat.gammaln(dof / 2));
|
2189 | }
|
2190 |
|
2191 |
|
2192 | return dof / x *
|
2193 | (jStat.noncentralt.cdf(x * Math.sqrt(1 + 2 / dof), dof+2, ncp) -
|
2194 | jStat.noncentralt.cdf(x, dof, ncp));
|
2195 | },
|
2196 |
|
2197 | cdf: function cdf(x, dof, ncp) {
|
2198 | var tol = 1e-14;
|
2199 | var min_iterations = 200;
|
2200 |
|
2201 | if (Math.abs(ncp) < tol)
|
2202 | return jStat.studentt.cdf(x, dof);
|
2203 |
|
2204 |
|
2205 | var flip = false;
|
2206 | if (x < 0) {
|
2207 | flip = true;
|
2208 | ncp = -ncp;
|
2209 | }
|
2210 |
|
2211 | var prob = jStat.normal.cdf(-ncp, 0, 1);
|
2212 | var value = tol + 1;
|
2213 |
|
2214 | var lastvalue = value;
|
2215 | var y = x * x / (x * x + dof);
|
2216 | var j = 0;
|
2217 | var p = Math.exp(-ncp * ncp / 2);
|
2218 | var q = Math.exp(-ncp * ncp / 2 - 0.5 * Math.log(2) -
|
2219 | jStat.gammaln(3 / 2)) * ncp;
|
2220 | while (j < min_iterations || lastvalue > tol || value > tol) {
|
2221 | lastvalue = value;
|
2222 | if (j > 0) {
|
2223 | p *= (ncp * ncp) / (2 * j);
|
2224 | q *= (ncp * ncp) / (2 * (j + 1 / 2));
|
2225 | }
|
2226 | value = p * jStat.beta.cdf(y, j + 0.5, dof / 2) +
|
2227 | q * jStat.beta.cdf(y, j+1, dof/2);
|
2228 | prob += 0.5 * value;
|
2229 | j++;
|
2230 | }
|
2231 |
|
2232 | return flip ? (1 - prob) : prob;
|
2233 | }
|
2234 | });
|
2235 |
|
2236 |
|
2237 |
|
2238 | jStat.extend(jStat.normal, {
|
2239 | pdf: function pdf(x, mean, std) {
|
2240 | return Math.exp(-0.5 * Math.log(2 * Math.PI) -
|
2241 | Math.log(std) - Math.pow(x - mean, 2) / (2 * std * std));
|
2242 | },
|
2243 |
|
2244 | cdf: function cdf(x, mean, std) {
|
2245 | return 0.5 * (1 + jStat.erf((x - mean) / Math.sqrt(2 * std * std)));
|
2246 | },
|
2247 |
|
2248 | inv: function(p, mean, std) {
|
2249 | return -1.41421356237309505 * std * jStat.erfcinv(2 * p) + mean;
|
2250 | },
|
2251 |
|
2252 | mean : function(mean/*, std*/) {
|
2253 | return mean;
|
2254 | },
|
2255 |
|
2256 | median: function median(mean/*, std*/) {
|
2257 | return mean;
|
2258 | },
|
2259 |
|
2260 | mode: function (mean/*, std*/) {
|
2261 | return mean;
|
2262 | },
|
2263 |
|
2264 | sample: function sample(mean, std) {
|
2265 | return jStat.randn() * std + mean;
|
2266 | },
|
2267 |
|
2268 | variance : function(mean, std) {
|
2269 | return std * std;
|
2270 | }
|
2271 | });
|
2272 |
|
2273 |
|
2274 |
|
2275 |
|
2276 | jStat.extend(jStat.pareto, {
|
2277 | pdf: function pdf(x, scale, shape) {
|
2278 | if (x < scale)
|
2279 | return 0;
|
2280 | return (shape * Math.pow(scale, shape)) / Math.pow(x, shape + 1);
|
2281 | },
|
2282 |
|
2283 | cdf: function cdf(x, scale, shape) {
|
2284 | if (x < scale)
|
2285 | return 0;
|
2286 | return 1 - Math.pow(scale / x, shape);
|
2287 | },
|
2288 |
|
2289 | inv: function inv(p, scale, shape) {
|
2290 | return scale / Math.pow(1 - p, 1 / shape);
|
2291 | },
|
2292 |
|
2293 | mean: function mean(scale, shape) {
|
2294 | if (shape <= 1)
|
2295 | return undefined;
|
2296 | return (shape * Math.pow(scale, shape)) / (shape - 1);
|
2297 | },
|
2298 |
|
2299 | median: function median(scale, shape) {
|
2300 | return scale * (shape * Math.SQRT2);
|
2301 | },
|
2302 |
|
2303 | mode: function mode(scale/*, shape*/) {
|
2304 | return scale;
|
2305 | },
|
2306 |
|
2307 | variance : function(scale, shape) {
|
2308 | if (shape <= 2)
|
2309 | return undefined;
|
2310 | return (scale*scale * shape) / (Math.pow(shape - 1, 2) * (shape - 2));
|
2311 | }
|
2312 | });
|
2313 |
|
2314 |
|
2315 |
|
2316 |
|
2317 | jStat.extend(jStat.studentt, {
|
2318 | pdf: function pdf(x, dof) {
|
2319 | dof = dof > 1e100 ? 1e100 : dof;
|
2320 | return (1/(Math.sqrt(dof) * jStat.betafn(0.5, dof/2))) *
|
2321 | Math.pow(1 + ((x * x) / dof), -((dof + 1) / 2));
|
2322 | },
|
2323 |
|
2324 | cdf: function cdf(x, dof) {
|
2325 | var dof2 = dof / 2;
|
2326 | return jStat.ibeta((x + Math.sqrt(x * x + dof)) /
|
2327 | (2 * Math.sqrt(x * x + dof)), dof2, dof2);
|
2328 | },
|
2329 |
|
2330 | inv: function(p, dof) {
|
2331 | var x = jStat.ibetainv(2 * Math.min(p, 1 - p), 0.5 * dof, 0.5);
|
2332 | x = Math.sqrt(dof * (1 - x) / x);
|
2333 | return (p > 0.5) ? x : -x;
|
2334 | },
|
2335 |
|
2336 | mean: function mean(dof) {
|
2337 | return (dof > 1) ? 0 : undefined;
|
2338 | },
|
2339 |
|
2340 | median: function median(/*dof*/) {
|
2341 | return 0;
|
2342 | },
|
2343 |
|
2344 | mode: function mode(/*dof*/) {
|
2345 | return 0;
|
2346 | },
|
2347 |
|
2348 | sample: function sample(dof) {
|
2349 | return jStat.randn() * Math.sqrt(dof / (2 * jStat.randg(dof / 2)));
|
2350 | },
|
2351 |
|
2352 | variance: function variance(dof) {
|
2353 | return (dof > 2) ? dof / (dof - 2) : (dof > 1) ? Infinity : undefined;
|
2354 | }
|
2355 | });
|
2356 |
|
2357 |
|
2358 |
|
2359 |
|
2360 | jStat.extend(jStat.weibull, {
|
2361 | pdf: function pdf(x, scale, shape) {
|
2362 | if (x < 0 || scale < 0 || shape < 0)
|
2363 | return 0;
|
2364 | return (shape / scale) * Math.pow((x / scale), (shape - 1)) *
|
2365 | Math.exp(-(Math.pow((x / scale), shape)));
|
2366 | },
|
2367 |
|
2368 | cdf: function cdf(x, scale, shape) {
|
2369 | return x < 0 ? 0 : 1 - Math.exp(-Math.pow((x / scale), shape));
|
2370 | },
|
2371 |
|
2372 | inv: function(p, scale, shape) {
|
2373 | return scale * Math.pow(-Math.log(1 - p), 1 / shape);
|
2374 | },
|
2375 |
|
2376 | mean : function(scale, shape) {
|
2377 | return scale * jStat.gammafn(1 + 1 / shape);
|
2378 | },
|
2379 |
|
2380 | median: function median(scale, shape) {
|
2381 | return scale * Math.pow(Math.log(2), 1 / shape);
|
2382 | },
|
2383 |
|
2384 | mode: function mode(scale, shape) {
|
2385 | if (shape <= 1)
|
2386 | return 0;
|
2387 | return scale * Math.pow((shape - 1) / shape, 1 / shape);
|
2388 | },
|
2389 |
|
2390 | sample: function sample(scale, shape) {
|
2391 | return scale * Math.pow(-Math.log(jStat._random_fn()), 1 / shape);
|
2392 | },
|
2393 |
|
2394 | variance: function variance(scale, shape) {
|
2395 | return scale * scale * jStat.gammafn(1 + 2 / shape) -
|
2396 | Math.pow(jStat.weibull.mean(scale, shape), 2);
|
2397 | }
|
2398 | });
|
2399 |
|
2400 |
|
2401 |
|
2402 |
|
2403 | jStat.extend(jStat.uniform, {
|
2404 | pdf: function pdf(x, a, b) {
|
2405 | return (x < a || x > b) ? 0 : 1 / (b - a);
|
2406 | },
|
2407 |
|
2408 | cdf: function cdf(x, a, b) {
|
2409 | if (x < a)
|
2410 | return 0;
|
2411 | else if (x < b)
|
2412 | return (x - a) / (b - a);
|
2413 | return 1;
|
2414 | },
|
2415 |
|
2416 | inv: function(p, a, b) {
|
2417 | return a + (p * (b - a));
|
2418 | },
|
2419 |
|
2420 | mean: function mean(a, b) {
|
2421 | return 0.5 * (a + b);
|
2422 | },
|
2423 |
|
2424 | median: function median(a, b) {
|
2425 | return jStat.mean(a, b);
|
2426 | },
|
2427 |
|
2428 | mode: function mode(/*a, b*/) {
|
2429 | throw new Error('mode is not yet implemented');
|
2430 | },
|
2431 |
|
2432 | sample: function sample(a, b) {
|
2433 | return (a / 2 + b / 2) + (b / 2 - a / 2) * (2 * jStat._random_fn() - 1);
|
2434 | },
|
2435 |
|
2436 | variance: function variance(a, b) {
|
2437 | return Math.pow(b - a, 2) / 12;
|
2438 | }
|
2439 | });
|
2440 |
|
2441 |
|
2442 |
|
2443 | function betinc(x, a, b, eps) {
|
2444 | var a0 = 0;
|
2445 | var b0 = 1;
|
2446 | var a1 = 1;
|
2447 | var b1 = 1;
|
2448 | var m9 = 0;
|
2449 | var a2 = 0;
|
2450 | var c9;
|
2451 |
|
2452 | while (Math.abs((a1 - a2) / a1) > eps) {
|
2453 | a2 = a1;
|
2454 | c9 = -(a + m9) * (a + b + m9) * x / (a + 2 * m9) / (a + 2 * m9 + 1);
|
2455 | a0 = a1 + c9 * a0;
|
2456 | b0 = b1 + c9 * b0;
|
2457 | m9 = m9 + 1;
|
2458 | c9 = m9 * (b - m9) * x / (a + 2 * m9 - 1) / (a + 2 * m9);
|
2459 | a1 = a0 + c9 * a1;
|
2460 | b1 = b0 + c9 * b1;
|
2461 | a0 = a0 / b1;
|
2462 | b0 = b0 / b1;
|
2463 | a1 = a1 / b1;
|
2464 | b1 = 1;
|
2465 | }
|
2466 |
|
2467 | return a1 / a;
|
2468 | }
|
2469 |
|
2470 |
|
2471 |
|
2472 | jStat.extend(jStat.binomial, {
|
2473 | pdf: function pdf(k, n, p) {
|
2474 | return (p === 0 || p === 1) ?
|
2475 | ((n * p) === k ? 1 : 0) :
|
2476 | jStat.combination(n, k) * Math.pow(p, k) * Math.pow(1 - p, n - k);
|
2477 | },
|
2478 |
|
2479 | cdf: function cdf(x, n, p) {
|
2480 | var betacdf;
|
2481 | var eps = 1e-10;
|
2482 |
|
2483 | if (x < 0)
|
2484 | return 0;
|
2485 | if (x >= n)
|
2486 | return 1;
|
2487 | if (p < 0 || p > 1 || n <= 0)
|
2488 | return NaN;
|
2489 |
|
2490 | x = Math.floor(x);
|
2491 | var z = p;
|
2492 | var a = x + 1;
|
2493 | var b = n - x;
|
2494 | var s = a + b;
|
2495 | var bt = Math.exp(jStat.gammaln(s) - jStat.gammaln(b) -
|
2496 | jStat.gammaln(a) + a * Math.log(z) + b * Math.log(1 - z));
|
2497 |
|
2498 | if (z < (a + 1) / (s + 2))
|
2499 | betacdf = bt * betinc(z, a, b, eps);
|
2500 | else
|
2501 | betacdf = 1 - bt * betinc(1 - z, b, a, eps);
|
2502 |
|
2503 | return Math.round((1 - betacdf) * (1 / eps)) / (1 / eps);
|
2504 | }
|
2505 | });
|
2506 |
|
2507 |
|
2508 |
|
2509 |
|
2510 | jStat.extend(jStat.negbin, {
|
2511 | pdf: function pdf(k, r, p) {
|
2512 | if (k !== k >>> 0)
|
2513 | return false;
|
2514 | if (k < 0)
|
2515 | return 0;
|
2516 | return jStat.combination(k + r - 1, r - 1) *
|
2517 | Math.pow(1 - p, k) * Math.pow(p, r);
|
2518 | },
|
2519 |
|
2520 | cdf: function cdf(x, r, p) {
|
2521 | var sum = 0,
|
2522 | k = 0;
|
2523 | if (x < 0) return 0;
|
2524 | for (; k <= x; k++) {
|
2525 | sum += jStat.negbin.pdf(k, r, p);
|
2526 | }
|
2527 | return sum;
|
2528 | }
|
2529 | });
|
2530 |
|
2531 |
|
2532 |
|
2533 |
|
2534 | jStat.extend(jStat.hypgeom, {
|
2535 | pdf: function pdf(k, N, m, n) {
|
2536 |
|
2537 |
|
2538 |
|
2539 |
|
2540 |
|
2541 |
|
2542 |
|
2543 |
|
2544 |
|
2545 | if(k !== k | 0) {
|
2546 | return false;
|
2547 | } else if(k < 0 || k < m - (N - n)) {
|
2548 |
|
2549 | return 0;
|
2550 | } else if(k > n || k > m) {
|
2551 |
|
2552 | return 0;
|
2553 | } else if (m * 2 > N) {
|
2554 |
|
2555 |
|
2556 | if(n * 2 > N) {
|
2557 |
|
2558 |
|
2559 | return jStat.hypgeom.pdf(N - m - n + k, N, N - m, N - n)
|
2560 | } else {
|
2561 |
|
2562 |
|
2563 | return jStat.hypgeom.pdf(n - k, N, N - m, n);
|
2564 | }
|
2565 |
|
2566 | } else if(n * 2 > N) {
|
2567 |
|
2568 |
|
2569 | return jStat.hypgeom.pdf(m - k, N, m, N - n);
|
2570 |
|
2571 | } else if(m < n) {
|
2572 |
|
2573 |
|
2574 | return jStat.hypgeom.pdf(k, N, n, m);
|
2575 | } else {
|
2576 |
|
2577 |
|
2578 |
|
2579 |
|
2580 |
|
2581 |
|
2582 |
|
2583 |
|
2584 |
|
2585 |
|
2586 |
|
2587 | var scaledPDF = 1;
|
2588 |
|
2589 |
|
2590 | var samplesDone = 0;
|
2591 |
|
2592 | for(var i = 0; i < k; i++) {
|
2593 |
|
2594 |
|
2595 | while(scaledPDF > 1 && samplesDone < n) {
|
2596 |
|
2597 |
|
2598 |
|
2599 | scaledPDF *= 1 - (m / (N - samplesDone));
|
2600 |
|
2601 |
|
2602 | samplesDone++;
|
2603 | }
|
2604 |
|
2605 |
|
2606 |
|
2607 | scaledPDF *= (n - i) * (m - i) / ((i + 1) * (N - m - n + i + 1));
|
2608 | }
|
2609 |
|
2610 | for(; samplesDone < n; samplesDone++) {
|
2611 |
|
2612 | scaledPDF *= 1 - (m / (N - samplesDone));
|
2613 | }
|
2614 |
|
2615 |
|
2616 | return Math.min(1, Math.max(0, scaledPDF));
|
2617 | }
|
2618 | },
|
2619 |
|
2620 | cdf: function cdf(x, N, m, n) {
|
2621 |
|
2622 |
|
2623 |
|
2624 |
|
2625 |
|
2626 |
|
2627 |
|
2628 |
|
2629 |
|
2630 |
|
2631 |
|
2632 | if(x < 0 || x < m - (N - n)) {
|
2633 |
|
2634 | return 0;
|
2635 | } else if(x >= n || x >= m) {
|
2636 |
|
2637 | return 1;
|
2638 | } else if (m * 2 > N) {
|
2639 |
|
2640 |
|
2641 | if(n * 2 > N) {
|
2642 |
|
2643 |
|
2644 | return jStat.hypgeom.cdf(N - m - n + x, N, N - m, N - n)
|
2645 | } else {
|
2646 |
|
2647 |
|
2648 | return 1 - jStat.hypgeom.cdf(n - x - 1, N, N - m, n);
|
2649 | }
|
2650 |
|
2651 | } else if(n * 2 > N) {
|
2652 |
|
2653 |
|
2654 | return 1 - jStat.hypgeom.cdf(m - x - 1, N, m, N - n);
|
2655 |
|
2656 | } else if(m < n) {
|
2657 |
|
2658 |
|
2659 | return jStat.hypgeom.cdf(x, N, n, m);
|
2660 | } else {
|
2661 |
|
2662 |
|
2663 |
|
2664 |
|
2665 |
|
2666 |
|
2667 |
|
2668 |
|
2669 |
|
2670 |
|
2671 | var scaledCDF = 1;
|
2672 |
|
2673 |
|
2674 |
|
2675 | var scaledPDF = 1;
|
2676 |
|
2677 |
|
2678 | var samplesDone = 0;
|
2679 |
|
2680 | for(var i = 0; i < x; i++) {
|
2681 |
|
2682 |
|
2683 | while(scaledCDF > 1 && samplesDone < n) {
|
2684 |
|
2685 |
|
2686 |
|
2687 | var factor = 1 - (m / (N - samplesDone));
|
2688 |
|
2689 | scaledPDF *= factor;
|
2690 | scaledCDF *= factor;
|
2691 |
|
2692 |
|
2693 | samplesDone++;
|
2694 | }
|
2695 |
|
2696 |
|
2697 |
|
2698 | scaledPDF *= (n - i) * (m - i) / ((i + 1) * (N - m - n + i + 1));
|
2699 |
|
2700 |
|
2701 | scaledCDF += scaledPDF;
|
2702 | }
|
2703 |
|
2704 | for(; samplesDone < n; samplesDone++) {
|
2705 |
|
2706 | scaledCDF *= 1 - (m / (N - samplesDone));
|
2707 | }
|
2708 |
|
2709 |
|
2710 | return Math.min(1, Math.max(0, scaledCDF));
|
2711 | }
|
2712 | }
|
2713 | });
|
2714 |
|
2715 |
|
2716 |
|
2717 |
|
2718 | jStat.extend(jStat.poisson, {
|
2719 | pdf: function pdf(k, l) {
|
2720 | if (l < 0 || (k % 1) !== 0 || k < 0) {
|
2721 | return 0;
|
2722 | }
|
2723 |
|
2724 | return Math.pow(l, k) * Math.exp(-l) / jStat.factorial(k);
|
2725 | },
|
2726 |
|
2727 | cdf: function cdf(x, l) {
|
2728 | var sumarr = [],
|
2729 | k = 0;
|
2730 | if (x < 0) return 0;
|
2731 | for (; k <= x; k++) {
|
2732 | sumarr.push(jStat.poisson.pdf(k, l));
|
2733 | }
|
2734 | return jStat.sum(sumarr);
|
2735 | },
|
2736 |
|
2737 | mean : function(l) {
|
2738 | return l;
|
2739 | },
|
2740 |
|
2741 | variance : function(l) {
|
2742 | return l;
|
2743 | },
|
2744 |
|
2745 | sample: function sample(l) {
|
2746 | var p = 1, k = 0, L = Math.exp(-l);
|
2747 | do {
|
2748 | k++;
|
2749 | p *= jStat._random_fn();
|
2750 | } while (p > L);
|
2751 | return k - 1;
|
2752 | }
|
2753 | });
|
2754 |
|
2755 |
|
2756 | jStat.extend(jStat.triangular, {
|
2757 | pdf: function pdf(x, a, b, c) {
|
2758 | if (b <= a || c < a || c > b) {
|
2759 | return NaN;
|
2760 | } else {
|
2761 | if (x < a || x > b) {
|
2762 | return 0;
|
2763 | } else if (x < c) {
|
2764 | return (2 * (x - a)) / ((b - a) * (c - a));
|
2765 | } else if (x === c) {
|
2766 | return (2 / (b - a));
|
2767 | } else {
|
2768 | return (2 * (b - x)) / ((b - a) * (b - c));
|
2769 | }
|
2770 | }
|
2771 | },
|
2772 |
|
2773 | cdf: function cdf(x, a, b, c) {
|
2774 | if (b <= a || c < a || c > b)
|
2775 | return NaN;
|
2776 | if (x <= a)
|
2777 | return 0;
|
2778 | else if (x >= b)
|
2779 | return 1;
|
2780 | if (x <= c)
|
2781 | return Math.pow(x - a, 2) / ((b - a) * (c - a));
|
2782 | else
|
2783 | return 1 - Math.pow(b - x, 2) / ((b - a) * (b - c));
|
2784 | },
|
2785 |
|
2786 | inv: function inv(p, a, b, c) {
|
2787 | if (b <= a || c < a || c > b) {
|
2788 | return NaN;
|
2789 | } else {
|
2790 | if (p <= ((c - a) / (b - a))) {
|
2791 | return a + (b - a) * Math.sqrt(p * ((c - a) / (b - a)));
|
2792 | } else {
|
2793 | return a + (b - a) * (1 - Math.sqrt((1 - p) * (1 - ((c - a) / (b - a)))));
|
2794 | }
|
2795 | }
|
2796 | },
|
2797 |
|
2798 | mean: function mean(a, b, c) {
|
2799 | return (a + b + c) / 3;
|
2800 | },
|
2801 |
|
2802 | median: function median(a, b, c) {
|
2803 | if (c <= (a + b) / 2) {
|
2804 | return b - Math.sqrt((b - a) * (b - c)) / Math.sqrt(2);
|
2805 | } else if (c > (a + b) / 2) {
|
2806 | return a + Math.sqrt((b - a) * (c - a)) / Math.sqrt(2);
|
2807 | }
|
2808 | },
|
2809 |
|
2810 | mode: function mode(a, b, c) {
|
2811 | return c;
|
2812 | },
|
2813 |
|
2814 | sample: function sample(a, b, c) {
|
2815 | var u = jStat._random_fn();
|
2816 | if (u < ((c - a) / (b - a)))
|
2817 | return a + Math.sqrt(u * (b - a) * (c - a))
|
2818 | return b - Math.sqrt((1 - u) * (b - a) * (b - c));
|
2819 | },
|
2820 |
|
2821 | variance: function variance(a, b, c) {
|
2822 | return (a * a + b * b + c * c - a * b - a * c - b * c) / 18;
|
2823 | }
|
2824 | });
|
2825 |
|
2826 |
|
2827 |
|
2828 | jStat.extend(jStat.arcsine, {
|
2829 | pdf: function pdf(x, a, b) {
|
2830 | if (b <= a) return NaN;
|
2831 |
|
2832 | return (x <= a || x >= b) ? 0 :
|
2833 | (2 / Math.PI) *
|
2834 | Math.pow(Math.pow(b - a, 2) -
|
2835 | Math.pow(2 * x - a - b, 2), -0.5);
|
2836 | },
|
2837 |
|
2838 | cdf: function cdf(x, a, b) {
|
2839 | if (x < a)
|
2840 | return 0;
|
2841 | else if (x < b)
|
2842 | return (2 / Math.PI) * Math.asin(Math.sqrt((x - a)/(b - a)));
|
2843 | return 1;
|
2844 | },
|
2845 |
|
2846 | inv: function(p, a, b) {
|
2847 | return a + (0.5 - 0.5 * Math.cos(Math.PI * p)) * (b - a);
|
2848 | },
|
2849 |
|
2850 | mean: function mean(a, b) {
|
2851 | if (b <= a) return NaN;
|
2852 | return (a + b) / 2;
|
2853 | },
|
2854 |
|
2855 | median: function median(a, b) {
|
2856 | if (b <= a) return NaN;
|
2857 | return (a + b) / 2;
|
2858 | },
|
2859 |
|
2860 | mode: function mode(/*a, b*/) {
|
2861 | throw new Error('mode is not yet implemented');
|
2862 | },
|
2863 |
|
2864 | sample: function sample(a, b) {
|
2865 | return ((a + b) / 2) + ((b - a) / 2) *
|
2866 | Math.sin(2 * Math.PI * jStat.uniform.sample(0, 1));
|
2867 | },
|
2868 |
|
2869 | variance: function variance(a, b) {
|
2870 | if (b <= a) return NaN;
|
2871 | return Math.pow(b - a, 2) / 8;
|
2872 | }
|
2873 | });
|
2874 |
|
2875 |
|
2876 | function laplaceSign(x) { return x / Math.abs(x); }
|
2877 |
|
2878 | jStat.extend(jStat.laplace, {
|
2879 | pdf: function pdf(x, mu, b) {
|
2880 | return (b <= 0) ? 0 : (Math.exp(-Math.abs(x - mu) / b)) / (2 * b);
|
2881 | },
|
2882 |
|
2883 | cdf: function cdf(x, mu, b) {
|
2884 | if (b <= 0) { return 0; }
|
2885 |
|
2886 | if(x < mu) {
|
2887 | return 0.5 * Math.exp((x - mu) / b);
|
2888 | } else {
|
2889 | return 1 - 0.5 * Math.exp(- (x - mu) / b);
|
2890 | }
|
2891 | },
|
2892 |
|
2893 | mean: function(mu/*, b*/) {
|
2894 | return mu;
|
2895 | },
|
2896 |
|
2897 | median: function(mu/*, b*/) {
|
2898 | return mu;
|
2899 | },
|
2900 |
|
2901 | mode: function(mu/*, b*/) {
|
2902 | return mu;
|
2903 | },
|
2904 |
|
2905 | variance: function(mu, b) {
|
2906 | return 2 * b * b;
|
2907 | },
|
2908 |
|
2909 | sample: function sample(mu, b) {
|
2910 | var u = jStat._random_fn() - 0.5;
|
2911 |
|
2912 | return mu - (b * laplaceSign(u) * Math.log(1 - (2 * Math.abs(u))));
|
2913 | }
|
2914 | });
|
2915 |
|
2916 | function tukeyWprob(w, rr, cc) {
|
2917 | var nleg = 12;
|
2918 | var ihalf = 6;
|
2919 |
|
2920 | var C1 = -30;
|
2921 | var C2 = -50;
|
2922 | var C3 = 60;
|
2923 | var bb = 8;
|
2924 | var wlar = 3;
|
2925 | var wincr1 = 2;
|
2926 | var wincr2 = 3;
|
2927 | var xleg = [
|
2928 | 0.981560634246719250690549090149,
|
2929 | 0.904117256370474856678465866119,
|
2930 | 0.769902674194304687036893833213,
|
2931 | 0.587317954286617447296702418941,
|
2932 | 0.367831498998180193752691536644,
|
2933 | 0.125233408511468915472441369464
|
2934 | ];
|
2935 | var aleg = [
|
2936 | 0.047175336386511827194615961485,
|
2937 | 0.106939325995318430960254718194,
|
2938 | 0.160078328543346226334652529543,
|
2939 | 0.203167426723065921749064455810,
|
2940 | 0.233492536538354808760849898925,
|
2941 | 0.249147045813402785000562436043
|
2942 | ];
|
2943 |
|
2944 | var qsqz = w * 0.5;
|
2945 |
|
2946 |
|
2947 |
|
2948 |
|
2949 | if (qsqz >= bb)
|
2950 | return 1.0;
|
2951 |
|
2952 |
|
2953 |
|
2954 |
|
2955 | var pr_w = 2 * jStat.normal.cdf(qsqz, 0, 1, 1, 0) - 1;
|
2956 |
|
2957 | if (pr_w >= Math.exp(C2 / cc))
|
2958 | pr_w = Math.pow(pr_w, cc);
|
2959 | else
|
2960 | pr_w = 0.0;
|
2961 |
|
2962 |
|
2963 |
|
2964 |
|
2965 | var wincr;
|
2966 | if (w > wlar)
|
2967 | wincr = wincr1;
|
2968 | else
|
2969 | wincr = wincr2;
|
2970 |
|
2971 |
|
2972 |
|
2973 |
|
2974 |
|
2975 |
|
2976 |
|
2977 |
|
2978 |
|
2979 | var blb = qsqz;
|
2980 | var binc = (bb - qsqz) / wincr;
|
2981 | var bub = blb + binc;
|
2982 | var einsum = 0.0;
|
2983 |
|
2984 |
|
2985 |
|
2986 | var cc1 = cc - 1.0;
|
2987 | for (var wi = 1; wi <= wincr; wi++) {
|
2988 | var elsum = 0.0;
|
2989 | var a = 0.5 * (bub + blb);
|
2990 |
|
2991 |
|
2992 |
|
2993 | var b = 0.5 * (bub - blb);
|
2994 |
|
2995 | for (var jj = 1; jj <= nleg; jj++) {
|
2996 | var j, xx;
|
2997 | if (ihalf < jj) {
|
2998 | j = (nleg - jj) + 1;
|
2999 | xx = xleg[j-1];
|
3000 | } else {
|
3001 | j = jj;
|
3002 | xx = -xleg[j-1];
|
3003 | }
|
3004 | var c = b * xx;
|
3005 | var ac = a + c;
|
3006 |
|
3007 |
|
3008 |
|
3009 |
|
3010 | var qexpo = ac * ac;
|
3011 | if (qexpo > C3)
|
3012 | break;
|
3013 |
|
3014 | var pplus = 2 * jStat.normal.cdf(ac, 0, 1, 1, 0);
|
3015 | var pminus= 2 * jStat.normal.cdf(ac, w, 1, 1, 0);
|
3016 |
|
3017 |
|
3018 |
|
3019 |
|
3020 | var rinsum = (pplus * 0.5) - (pminus * 0.5);
|
3021 | if (rinsum >= Math.exp(C1 / cc1)) {
|
3022 | rinsum = (aleg[j-1] * Math.exp(-(0.5 * qexpo))) * Math.pow(rinsum, cc1);
|
3023 | elsum += rinsum;
|
3024 | }
|
3025 | }
|
3026 | elsum *= (((2.0 * b) * cc) / Math.sqrt(2 * Math.PI));
|
3027 | einsum += elsum;
|
3028 | blb = bub;
|
3029 | bub += binc;
|
3030 | }
|
3031 |
|
3032 |
|
3033 | pr_w += einsum;
|
3034 | if (pr_w <= Math.exp(C1 / rr))
|
3035 | return 0;
|
3036 |
|
3037 | pr_w = Math.pow(pr_w, rr);
|
3038 | if (pr_w >= 1)
|
3039 | return 1;
|
3040 | return pr_w;
|
3041 | }
|
3042 |
|
3043 | function tukeyQinv(p, c, v) {
|
3044 | var p0 = 0.322232421088;
|
3045 | var q0 = 0.993484626060e-01;
|
3046 | var p1 = -1.0;
|
3047 | var q1 = 0.588581570495;
|
3048 | var p2 = -0.342242088547;
|
3049 | var q2 = 0.531103462366;
|
3050 | var p3 = -0.204231210125;
|
3051 | var q3 = 0.103537752850;
|
3052 | var p4 = -0.453642210148e-04;
|
3053 | var q4 = 0.38560700634e-02;
|
3054 | var c1 = 0.8832;
|
3055 | var c2 = 0.2368;
|
3056 | var c3 = 1.214;
|
3057 | var c4 = 1.208;
|
3058 | var c5 = 1.4142;
|
3059 | var vmax = 120.0;
|
3060 |
|
3061 | var ps = 0.5 - 0.5 * p;
|
3062 | var yi = Math.sqrt(Math.log(1.0 / (ps * ps)));
|
3063 | var t = yi + (((( yi * p4 + p3) * yi + p2) * yi + p1) * yi + p0)
|
3064 | / (((( yi * q4 + q3) * yi + q2) * yi + q1) * yi + q0);
|
3065 | if (v < vmax) t += (t * t * t + t) / v / 4.0;
|
3066 | var q = c1 - c2 * t;
|
3067 | if (v < vmax) q += -c3 / v + c4 * t / v;
|
3068 | return t * (q * Math.log(c - 1.0) + c5);
|
3069 | }
|
3070 |
|
3071 | jStat.extend(jStat.tukey, {
|
3072 | cdf: function cdf(q, nmeans, df) {
|
3073 |
|
3074 | var rr = 1;
|
3075 | var cc = nmeans;
|
3076 |
|
3077 | var nlegq = 16;
|
3078 | var ihalfq = 8;
|
3079 |
|
3080 | var eps1 = -30.0;
|
3081 | var eps2 = 1.0e-14;
|
3082 | var dhaf = 100.0;
|
3083 | var dquar = 800.0;
|
3084 | var deigh = 5000.0;
|
3085 | var dlarg = 25000.0;
|
3086 | var ulen1 = 1.0;
|
3087 | var ulen2 = 0.5;
|
3088 | var ulen3 = 0.25;
|
3089 | var ulen4 = 0.125;
|
3090 | var xlegq = [
|
3091 | 0.989400934991649932596154173450,
|
3092 | 0.944575023073232576077988415535,
|
3093 | 0.865631202387831743880467897712,
|
3094 | 0.755404408355003033895101194847,
|
3095 | 0.617876244402643748446671764049,
|
3096 | 0.458016777657227386342419442984,
|
3097 | 0.281603550779258913230460501460,
|
3098 | 0.950125098376374401853193354250e-1
|
3099 | ];
|
3100 | var alegq = [
|
3101 | 0.271524594117540948517805724560e-1,
|
3102 | 0.622535239386478928628438369944e-1,
|
3103 | 0.951585116824927848099251076022e-1,
|
3104 | 0.124628971255533872052476282192,
|
3105 | 0.149595988816576732081501730547,
|
3106 | 0.169156519395002538189312079030,
|
3107 | 0.182603415044923588866763667969,
|
3108 | 0.189450610455068496285396723208
|
3109 | ];
|
3110 |
|
3111 | if (q <= 0)
|
3112 | return 0;
|
3113 |
|
3114 |
|
3115 |
|
3116 |
|
3117 | if (df < 2 || rr < 1 || cc < 2) return NaN;
|
3118 |
|
3119 | if (!Number.isFinite(q))
|
3120 | return 1;
|
3121 |
|
3122 | if (df > dlarg)
|
3123 | return tukeyWprob(q, rr, cc);
|
3124 |
|
3125 |
|
3126 |
|
3127 | var f2 = df * 0.5;
|
3128 | var f2lf = ((f2 * Math.log(df)) - (df * Math.log(2))) - jStat.gammaln(f2);
|
3129 | var f21 = f2 - 1.0;
|
3130 |
|
3131 |
|
3132 |
|
3133 |
|
3134 |
|
3135 | var ff4 = df * 0.25;
|
3136 | var ulen;
|
3137 | if (df <= dhaf) ulen = ulen1;
|
3138 | else if (df <= dquar) ulen = ulen2;
|
3139 | else if (df <= deigh) ulen = ulen3;
|
3140 | else ulen = ulen4;
|
3141 |
|
3142 | f2lf += Math.log(ulen);
|
3143 |
|
3144 |
|
3145 |
|
3146 | var ans = 0.0;
|
3147 |
|
3148 | for (var i = 1; i <= 50; i++) {
|
3149 | var otsum = 0.0;
|
3150 |
|
3151 |
|
3152 |
|
3153 |
|
3154 | var twa1 = (2 * i - 1) * ulen;
|
3155 |
|
3156 | for (var jj = 1; jj <= nlegq; jj++) {
|
3157 | var j, t1;
|
3158 | if (ihalfq < jj) {
|
3159 | j = jj - ihalfq - 1;
|
3160 | t1 = (f2lf + (f21 * Math.log(twa1 + (xlegq[j] * ulen))))
|
3161 | - (((xlegq[j] * ulen) + twa1) * ff4);
|
3162 | } else {
|
3163 | j = jj - 1;
|
3164 | t1 = (f2lf + (f21 * Math.log(twa1 - (xlegq[j] * ulen))))
|
3165 | + (((xlegq[j] * ulen) - twa1) * ff4);
|
3166 | }
|
3167 |
|
3168 |
|
3169 | var qsqz;
|
3170 | if (t1 >= eps1) {
|
3171 | if (ihalfq < jj) {
|
3172 | qsqz = q * Math.sqrt(((xlegq[j] * ulen) + twa1) * 0.5);
|
3173 | } else {
|
3174 | qsqz = q * Math.sqrt(((-(xlegq[j] * ulen)) + twa1) * 0.5);
|
3175 | }
|
3176 |
|
3177 |
|
3178 |
|
3179 | var wprb = tukeyWprob(qsqz, rr, cc);
|
3180 | var rotsum = (wprb * alegq[j]) * Math.exp(t1);
|
3181 | otsum += rotsum;
|
3182 | }
|
3183 |
|
3184 |
|
3185 | }
|
3186 |
|
3187 |
|
3188 |
|
3189 |
|
3190 | if (i * ulen >= 1.0 && otsum <= eps2)
|
3191 | break;
|
3192 |
|
3193 |
|
3194 |
|
3195 |
|
3196 | ans += otsum;
|
3197 | }
|
3198 |
|
3199 | if (otsum > eps2) {
|
3200 | throw new Error('tukey.cdf failed to converge');
|
3201 | }
|
3202 | if (ans > 1)
|
3203 | ans = 1;
|
3204 | return ans;
|
3205 | },
|
3206 |
|
3207 | inv: function(p, nmeans, df) {
|
3208 |
|
3209 | var rr = 1;
|
3210 | var cc = nmeans;
|
3211 |
|
3212 | var eps = 0.0001;
|
3213 | var maxiter = 50;
|
3214 |
|
3215 |
|
3216 | if (df < 2 || rr < 1 || cc < 2) return NaN;
|
3217 |
|
3218 | if (p < 0 || p > 1) return NaN;
|
3219 | if (p === 0) return 0;
|
3220 | if (p === 1) return Infinity;
|
3221 |
|
3222 |
|
3223 |
|
3224 | var x0 = tukeyQinv(p, cc, df);
|
3225 |
|
3226 |
|
3227 |
|
3228 | var valx0 = jStat.tukey.cdf(x0, nmeans, df) - p;
|
3229 |
|
3230 |
|
3231 |
|
3232 |
|
3233 |
|
3234 |
|
3235 | var x1;
|
3236 | if (valx0 > 0.0)
|
3237 | x1 = Math.max(0.0, x0 - 1.0);
|
3238 | else
|
3239 | x1 = x0 + 1.0;
|
3240 | var valx1 = jStat.tukey.cdf(x1, nmeans, df) - p;
|
3241 |
|
3242 |
|
3243 |
|
3244 | var ans;
|
3245 | for(var iter = 1; iter < maxiter; iter++) {
|
3246 | ans = x1 - ((valx1 * (x1 - x0)) / (valx1 - valx0));
|
3247 | valx0 = valx1;
|
3248 |
|
3249 |
|
3250 |
|
3251 | x0 = x1;
|
3252 | if (ans < 0.0) {
|
3253 | ans = 0.0;
|
3254 | valx1 = -p;
|
3255 | }
|
3256 |
|
3257 |
|
3258 | valx1 = jStat.tukey.cdf(ans, nmeans, df) - p;
|
3259 | x1 = ans;
|
3260 |
|
3261 |
|
3262 |
|
3263 |
|
3264 | var xabs = Math.abs(x1 - x0);
|
3265 | if (xabs < eps)
|
3266 | return ans;
|
3267 | }
|
3268 |
|
3269 | throw new Error('tukey.inv failed to converge');
|
3270 | }
|
3271 | });
|
3272 |
|
3273 | }(jStat, Math));
|
3274 |
|
3275 |
|
3276 |
|
3277 | (function(jStat, Math) {
|
3278 |
|
3279 | var push = Array.prototype.push;
|
3280 | var isArray = jStat.utils.isArray;
|
3281 |
|
3282 | function isUsable(arg) {
|
3283 | return isArray(arg) || arg instanceof jStat;
|
3284 | }
|
3285 |
|
3286 | jStat.extend({
|
3287 |
|
3288 |
|
3289 | add: function add(arr, arg) {
|
3290 |
|
3291 | if (isUsable(arg)) {
|
3292 | if (!isUsable(arg[0])) arg = [ arg ];
|
3293 | return jStat.map(arr, function(value, row, col) {
|
3294 | return value + arg[row][col];
|
3295 | });
|
3296 | }
|
3297 | return jStat.map(arr, function(value) { return value + arg; });
|
3298 | },
|
3299 |
|
3300 |
|
3301 | subtract: function subtract(arr, arg) {
|
3302 |
|
3303 | if (isUsable(arg)) {
|
3304 | if (!isUsable(arg[0])) arg = [ arg ];
|
3305 | return jStat.map(arr, function(value, row, col) {
|
3306 | return value - arg[row][col] || 0;
|
3307 | });
|
3308 | }
|
3309 | return jStat.map(arr, function(value) { return value - arg; });
|
3310 | },
|
3311 |
|
3312 |
|
3313 | divide: function divide(arr, arg) {
|
3314 | if (isUsable(arg)) {
|
3315 | if (!isUsable(arg[0])) arg = [ arg ];
|
3316 | return jStat.multiply(arr, jStat.inv(arg));
|
3317 | }
|
3318 | return jStat.map(arr, function(value) { return value / arg; });
|
3319 | },
|
3320 |
|
3321 |
|
3322 | multiply: function multiply(arr, arg) {
|
3323 | var row, col, nrescols, sum, nrow, ncol, res, rescols;
|
3324 |
|
3325 | if (arr.length === undefined && arg.length === undefined) {
|
3326 | return arr * arg;
|
3327 | }
|
3328 | nrow = arr.length,
|
3329 | ncol = arr[0].length,
|
3330 | res = jStat.zeros(nrow, nrescols = (isUsable(arg)) ? arg[0].length : ncol),
|
3331 | rescols = 0;
|
3332 | if (isUsable(arg)) {
|
3333 | for (; rescols < nrescols; rescols++) {
|
3334 | for (row = 0; row < nrow; row++) {
|
3335 | sum = 0;
|
3336 | for (col = 0; col < ncol; col++)
|
3337 | sum += arr[row][col] * arg[col][rescols];
|
3338 | res[row][rescols] = sum;
|
3339 | }
|
3340 | }
|
3341 | return (nrow === 1 && rescols === 1) ? res[0][0] : res;
|
3342 | }
|
3343 | return jStat.map(arr, function(value) { return value * arg; });
|
3344 | },
|
3345 |
|
3346 |
|
3347 |
|
3348 |
|
3349 |
|
3350 |
|
3351 | outer:function outer(A, B) {
|
3352 | return jStat.multiply(A.map(function(t){ return [t] }), [B]);
|
3353 | },
|
3354 |
|
3355 |
|
3356 |
|
3357 | dot: function dot(arr, arg) {
|
3358 | if (!isUsable(arr[0])) arr = [ arr ];
|
3359 | if (!isUsable(arg[0])) arg = [ arg ];
|
3360 |
|
3361 | var left = (arr[0].length === 1 && arr.length !== 1) ? jStat.transpose(arr) : arr,
|
3362 | right = (arg[0].length === 1 && arg.length !== 1) ? jStat.transpose(arg) : arg,
|
3363 | res = [],
|
3364 | row = 0,
|
3365 | nrow = left.length,
|
3366 | ncol = left[0].length,
|
3367 | sum, col;
|
3368 | for (; row < nrow; row++) {
|
3369 | res[row] = [];
|
3370 | sum = 0;
|
3371 | for (col = 0; col < ncol; col++)
|
3372 | sum += left[row][col] * right[row][col];
|
3373 | res[row] = sum;
|
3374 | }
|
3375 | return (res.length === 1) ? res[0] : res;
|
3376 | },
|
3377 |
|
3378 |
|
3379 | pow: function pow(arr, arg) {
|
3380 | return jStat.map(arr, function(value) { return Math.pow(value, arg); });
|
3381 | },
|
3382 |
|
3383 |
|
3384 | exp: function exp(arr) {
|
3385 | return jStat.map(arr, function(value) { return Math.exp(value); });
|
3386 | },
|
3387 |
|
3388 |
|
3389 | log: function exp(arr) {
|
3390 | return jStat.map(arr, function(value) { return Math.log(value); });
|
3391 | },
|
3392 |
|
3393 |
|
3394 | abs: function abs(arr) {
|
3395 | return jStat.map(arr, function(value) { return Math.abs(value); });
|
3396 | },
|
3397 |
|
3398 |
|
3399 |
|
3400 | norm: function norm(arr, p) {
|
3401 | var nnorm = 0,
|
3402 | i = 0;
|
3403 |
|
3404 | if (isNaN(p)) p = 2;
|
3405 |
|
3406 | if (isUsable(arr[0])) arr = arr[0];
|
3407 |
|
3408 | for (; i < arr.length; i++) {
|
3409 | nnorm += Math.pow(Math.abs(arr[i]), p);
|
3410 | }
|
3411 | return Math.pow(nnorm, 1 / p);
|
3412 | },
|
3413 |
|
3414 |
|
3415 |
|
3416 | angle: function angle(arr, arg) {
|
3417 | return Math.acos(jStat.dot(arr, arg) / (jStat.norm(arr) * jStat.norm(arg)));
|
3418 | },
|
3419 |
|
3420 |
|
3421 |
|
3422 | aug: function aug(a, b) {
|
3423 | var newarr = [];
|
3424 | var i;
|
3425 | for (i = 0; i < a.length; i++) {
|
3426 | newarr.push(a[i].slice());
|
3427 | }
|
3428 | for (i = 0; i < newarr.length; i++) {
|
3429 | push.apply(newarr[i], b[i]);
|
3430 | }
|
3431 | return newarr;
|
3432 | },
|
3433 |
|
3434 |
|
3435 |
|
3436 |
|
3437 | inv: function inv(a) {
|
3438 | var rows = a.length;
|
3439 | var cols = a[0].length;
|
3440 | var b = jStat.identity(rows, cols);
|
3441 | var c = jStat.gauss_jordan(a, b);
|
3442 | var result = [];
|
3443 | var i = 0;
|
3444 | var j;
|
3445 |
|
3446 |
|
3447 | for (; i < rows; i++) {
|
3448 | result[i] = [];
|
3449 | for (j = cols; j < c[0].length; j++)
|
3450 | result[i][j - cols] = c[i][j];
|
3451 | }
|
3452 | return result;
|
3453 | },
|
3454 |
|
3455 |
|
3456 | det: function det(a) {
|
3457 | var alen = a.length,
|
3458 | alend = alen * 2,
|
3459 | vals = new Array(alend),
|
3460 | rowshift = alen - 1,
|
3461 | colshift = alend - 1,
|
3462 | mrow = rowshift - alen + 1,
|
3463 | mcol = colshift,
|
3464 | i = 0,
|
3465 | result = 0,
|
3466 | j;
|
3467 |
|
3468 | if (alen === 2) {
|
3469 | return a[0][0] * a[1][1] - a[0][1] * a[1][0];
|
3470 | }
|
3471 | for (; i < alend; i++) {
|
3472 | vals[i] = 1;
|
3473 | }
|
3474 | for (i = 0; i < alen; i++) {
|
3475 | for (j = 0; j < alen; j++) {
|
3476 | vals[(mrow < 0) ? mrow + alen : mrow ] *= a[i][j];
|
3477 | vals[(mcol < alen) ? mcol + alen : mcol ] *= a[i][j];
|
3478 | mrow++;
|
3479 | mcol--;
|
3480 | }
|
3481 | mrow = --rowshift - alen + 1;
|
3482 | mcol = --colshift;
|
3483 | }
|
3484 | for (i = 0; i < alen; i++) {
|
3485 | result += vals[i];
|
3486 | }
|
3487 | for (; i < alend; i++) {
|
3488 | result -= vals[i];
|
3489 | }
|
3490 | return result;
|
3491 | },
|
3492 |
|
3493 | gauss_elimination: function gauss_elimination(a, b) {
|
3494 | var i = 0,
|
3495 | j = 0,
|
3496 | n = a.length,
|
3497 | m = a[0].length,
|
3498 | factor = 1,
|
3499 | sum = 0,
|
3500 | x = [],
|
3501 | maug, pivot, temp, k;
|
3502 | a = jStat.aug(a, b);
|
3503 | maug = a[0].length;
|
3504 | for(i = 0; i < n; i++) {
|
3505 | pivot = a[i][i];
|
3506 | j = i;
|
3507 | for (k = i + 1; k < m; k++) {
|
3508 | if (pivot < Math.abs(a[k][i])) {
|
3509 | pivot = a[k][i];
|
3510 | j = k;
|
3511 | }
|
3512 | }
|
3513 | if (j != i) {
|
3514 | for(k = 0; k < maug; k++) {
|
3515 | temp = a[i][k];
|
3516 | a[i][k] = a[j][k];
|
3517 | a[j][k] = temp;
|
3518 | }
|
3519 | }
|
3520 | for (j = i + 1; j < n; j++) {
|
3521 | factor = a[j][i] / a[i][i];
|
3522 | for(k = i; k < maug; k++) {
|
3523 | a[j][k] = a[j][k] - factor * a[i][k];
|
3524 | }
|
3525 | }
|
3526 | }
|
3527 | for (i = n - 1; i >= 0; i--) {
|
3528 | sum = 0;
|
3529 | for (j = i + 1; j<= n - 1; j++) {
|
3530 | sum = sum + x[j] * a[i][j];
|
3531 | }
|
3532 | x[i] =(a[i][maug - 1] - sum) / a[i][i];
|
3533 | }
|
3534 | return x;
|
3535 | },
|
3536 |
|
3537 | gauss_jordan: function gauss_jordan(a, b) {
|
3538 | var m = jStat.aug(a, b);
|
3539 | var h = m.length;
|
3540 | var w = m[0].length;
|
3541 | var c = 0;
|
3542 | var x, y, y2;
|
3543 |
|
3544 | for (y = 0; y < h; y++) {
|
3545 | var maxrow = y;
|
3546 | for (y2 = y+1; y2 < h; y2++) {
|
3547 | if (Math.abs(m[y2][y]) > Math.abs(m[maxrow][y]))
|
3548 | maxrow = y2;
|
3549 | }
|
3550 | var tmp = m[y];
|
3551 | m[y] = m[maxrow];
|
3552 | m[maxrow] = tmp
|
3553 | for (y2 = y+1; y2 < h; y2++) {
|
3554 | c = m[y2][y] / m[y][y];
|
3555 | for (x = y; x < w; x++) {
|
3556 | m[y2][x] -= m[y][x] * c;
|
3557 | }
|
3558 | }
|
3559 | }
|
3560 |
|
3561 | for (y = h-1; y >= 0; y--) {
|
3562 | c = m[y][y];
|
3563 | for (y2 = 0; y2 < y; y2++) {
|
3564 | for (x = w-1; x > y-1; x--) {
|
3565 | m[y2][x] -= m[y][x] * m[y2][y] / c;
|
3566 | }
|
3567 | }
|
3568 | m[y][y] /= c;
|
3569 | for (x = h; x < w; x++) {
|
3570 | m[y][x] /= c;
|
3571 | }
|
3572 | }
|
3573 | return m;
|
3574 | },
|
3575 |
|
3576 |
|
3577 |
|
3578 |
|
3579 |
|
3580 |
|
3581 |
|
3582 |
|
3583 |
|
3584 |
|
3585 |
|
3586 | triaUpSolve: function triaUpSolve(A, b) {
|
3587 | var size = A[0].length;
|
3588 | var x = jStat.zeros(1, size)[0];
|
3589 | var parts;
|
3590 | var matrix_mode = false;
|
3591 |
|
3592 | if (b[0].length != undefined) {
|
3593 | b = b.map(function(i){ return i[0] });
|
3594 | matrix_mode = true;
|
3595 | }
|
3596 |
|
3597 | jStat.arange(size - 1, -1, -1).forEach(function(i) {
|
3598 | parts = jStat.arange(i + 1, size).map(function(j) {
|
3599 | return x[j] * A[i][j];
|
3600 | });
|
3601 | x[i] = (b[i] - jStat.sum(parts)) / A[i][i];
|
3602 | });
|
3603 |
|
3604 | if (matrix_mode)
|
3605 | return x.map(function(i){ return [i] });
|
3606 | return x;
|
3607 | },
|
3608 |
|
3609 | triaLowSolve: function triaLowSolve(A, b) {
|
3610 |
|
3611 | var size = A[0].length;
|
3612 | var x = jStat.zeros(1, size)[0];
|
3613 | var parts;
|
3614 |
|
3615 | var matrix_mode=false;
|
3616 | if (b[0].length != undefined) {
|
3617 | b = b.map(function(i){ return i[0] });
|
3618 | matrix_mode = true;
|
3619 | }
|
3620 |
|
3621 | jStat.arange(size).forEach(function(i) {
|
3622 | parts = jStat.arange(i).map(function(j) {
|
3623 | return A[i][j] * x[j];
|
3624 | });
|
3625 | x[i] = (b[i] - jStat.sum(parts)) / A[i][i];
|
3626 | })
|
3627 |
|
3628 | if (matrix_mode)
|
3629 | return x.map(function(i){ return [i] });
|
3630 | return x;
|
3631 | },
|
3632 |
|
3633 |
|
3634 |
|
3635 |
|
3636 |
|
3637 |
|
3638 | lu: function lu(A) {
|
3639 | var size = A.length;
|
3640 |
|
3641 | var L = jStat.identity(size);
|
3642 | var R = jStat.zeros(A.length, A[0].length);
|
3643 | var parts;
|
3644 | jStat.arange(size).forEach(function(t) {
|
3645 | R[0][t] = A[0][t];
|
3646 | });
|
3647 | jStat.arange(1, size).forEach(function(l) {
|
3648 | jStat.arange(l).forEach(function(i) {
|
3649 | parts = jStat.arange(i).map(function(jj) {
|
3650 | return L[l][jj] * R[jj][i];
|
3651 | });
|
3652 | L[l][i] = (A[l][i] - jStat.sum(parts)) / R[i][i];
|
3653 | });
|
3654 | jStat.arange(l, size).forEach(function(j) {
|
3655 | parts = jStat.arange(l).map(function(jj) {
|
3656 | return L[l][jj] * R[jj][j];
|
3657 | });
|
3658 | R[l][j] = A[parts.length][j] - jStat.sum(parts);
|
3659 | });
|
3660 | });
|
3661 | return [L, R];
|
3662 | },
|
3663 |
|
3664 |
|
3665 |
|
3666 |
|
3667 | cholesky: function cholesky(A) {
|
3668 | var size = A.length;
|
3669 | var T = jStat.zeros(A.length, A[0].length);
|
3670 | var parts;
|
3671 | jStat.arange(size).forEach(function(i) {
|
3672 | parts = jStat.arange(i).map(function(t) {
|
3673 | return Math.pow(T[i][t],2);
|
3674 | });
|
3675 | T[i][i] = Math.sqrt(A[i][i] - jStat.sum(parts));
|
3676 | jStat.arange(i + 1, size).forEach(function(j) {
|
3677 | parts = jStat.arange(i).map(function(t) {
|
3678 | return T[i][t] * T[j][t];
|
3679 | });
|
3680 | T[j][i] = (A[i][j] - jStat.sum(parts)) / T[i][i];
|
3681 | });
|
3682 | });
|
3683 | return T;
|
3684 | },
|
3685 |
|
3686 |
|
3687 | gauss_jacobi: function gauss_jacobi(a, b, x, r) {
|
3688 | var i = 0;
|
3689 | var j = 0;
|
3690 | var n = a.length;
|
3691 | var l = [];
|
3692 | var u = [];
|
3693 | var d = [];
|
3694 | var xv, c, h, xk;
|
3695 | for (; i < n; i++) {
|
3696 | l[i] = [];
|
3697 | u[i] = [];
|
3698 | d[i] = [];
|
3699 | for (j = 0; j < n; j++) {
|
3700 | if (i > j) {
|
3701 | l[i][j] = a[i][j];
|
3702 | u[i][j] = d[i][j] = 0;
|
3703 | } else if (i < j) {
|
3704 | u[i][j] = a[i][j];
|
3705 | l[i][j] = d[i][j] = 0;
|
3706 | } else {
|
3707 | d[i][j] = a[i][j];
|
3708 | l[i][j] = u[i][j] = 0;
|
3709 | }
|
3710 | }
|
3711 | }
|
3712 | h = jStat.multiply(jStat.multiply(jStat.inv(d), jStat.add(l, u)), -1);
|
3713 | c = jStat.multiply(jStat.inv(d), b);
|
3714 | xv = x;
|
3715 | xk = jStat.add(jStat.multiply(h, x), c);
|
3716 | i = 2;
|
3717 | while (Math.abs(jStat.norm(jStat.subtract(xk,xv))) > r) {
|
3718 | xv = xk;
|
3719 | xk = jStat.add(jStat.multiply(h, xv), c);
|
3720 | i++;
|
3721 | }
|
3722 | return xk;
|
3723 | },
|
3724 |
|
3725 | gauss_seidel: function gauss_seidel(a, b, x, r) {
|
3726 | var i = 0;
|
3727 | var n = a.length;
|
3728 | var l = [];
|
3729 | var u = [];
|
3730 | var d = [];
|
3731 | var j, xv, c, h, xk;
|
3732 | for (; i < n; i++) {
|
3733 | l[i] = [];
|
3734 | u[i] = [];
|
3735 | d[i] = [];
|
3736 | for (j = 0; j < n; j++) {
|
3737 | if (i > j) {
|
3738 | l[i][j] = a[i][j];
|
3739 | u[i][j] = d[i][j] = 0;
|
3740 | } else if (i < j) {
|
3741 | u[i][j] = a[i][j];
|
3742 | l[i][j] = d[i][j] = 0;
|
3743 | } else {
|
3744 | d[i][j] = a[i][j];
|
3745 | l[i][j] = u[i][j] = 0;
|
3746 | }
|
3747 | }
|
3748 | }
|
3749 | h = jStat.multiply(jStat.multiply(jStat.inv(jStat.add(d, l)), u), -1);
|
3750 | c = jStat.multiply(jStat.inv(jStat.add(d, l)), b);
|
3751 | xv = x;
|
3752 | xk = jStat.add(jStat.multiply(h, x), c);
|
3753 | i = 2;
|
3754 | while (Math.abs(jStat.norm(jStat.subtract(xk, xv))) > r) {
|
3755 | xv = xk;
|
3756 | xk = jStat.add(jStat.multiply(h, xv), c);
|
3757 | i = i + 1;
|
3758 | }
|
3759 | return xk;
|
3760 | },
|
3761 |
|
3762 | SOR: function SOR(a, b, x, r, w) {
|
3763 | var i = 0;
|
3764 | var n = a.length;
|
3765 | var l = [];
|
3766 | var u = [];
|
3767 | var d = [];
|
3768 | var j, xv, c, h, xk;
|
3769 | for (; i < n; i++) {
|
3770 | l[i] = [];
|
3771 | u[i] = [];
|
3772 | d[i] = [];
|
3773 | for (j = 0; j < n; j++) {
|
3774 | if (i > j) {
|
3775 | l[i][j] = a[i][j];
|
3776 | u[i][j] = d[i][j] = 0;
|
3777 | } else if (i < j) {
|
3778 | u[i][j] = a[i][j];
|
3779 | l[i][j] = d[i][j] = 0;
|
3780 | } else {
|
3781 | d[i][j] = a[i][j];
|
3782 | l[i][j] = u[i][j] = 0;
|
3783 | }
|
3784 | }
|
3785 | }
|
3786 | h = jStat.multiply(jStat.inv(jStat.add(d, jStat.multiply(l, w))),
|
3787 | jStat.subtract(jStat.multiply(d, 1 - w),
|
3788 | jStat.multiply(u, w)));
|
3789 | c = jStat.multiply(jStat.multiply(jStat.inv(jStat.add(d,
|
3790 | jStat.multiply(l, w))), b), w);
|
3791 | xv = x;
|
3792 | xk = jStat.add(jStat.multiply(h, x), c);
|
3793 | i = 2;
|
3794 | while (Math.abs(jStat.norm(jStat.subtract(xk, xv))) > r) {
|
3795 | xv = xk;
|
3796 | xk = jStat.add(jStat.multiply(h, xv), c);
|
3797 | i++;
|
3798 | }
|
3799 | return xk;
|
3800 | },
|
3801 |
|
3802 | householder: function householder(a) {
|
3803 | var m = a.length;
|
3804 | var n = a[0].length;
|
3805 | var i = 0;
|
3806 | var w = [];
|
3807 | var p = [];
|
3808 | var alpha, r, k, j, factor;
|
3809 | for (; i < m - 1; i++) {
|
3810 | alpha = 0;
|
3811 | for (j = i + 1; j < n; j++)
|
3812 | alpha += (a[j][i] * a[j][i]);
|
3813 | factor = (a[i + 1][i] > 0) ? -1 : 1;
|
3814 | alpha = factor * Math.sqrt(alpha);
|
3815 | r = Math.sqrt((((alpha * alpha) - a[i + 1][i] * alpha) / 2));
|
3816 | w = jStat.zeros(m, 1);
|
3817 | w[i + 1][0] = (a[i + 1][i] - alpha) / (2 * r);
|
3818 | for (k = i + 2; k < m; k++) w[k][0] = a[k][i] / (2 * r);
|
3819 | p = jStat.subtract(jStat.identity(m, n),
|
3820 | jStat.multiply(jStat.multiply(w, jStat.transpose(w)), 2));
|
3821 | a = jStat.multiply(p, jStat.multiply(a, p));
|
3822 | }
|
3823 | return a;
|
3824 | },
|
3825 |
|
3826 |
|
3827 |
|
3828 |
|
3829 | QR: (function() {
|
3830 |
|
3831 |
|
3832 |
|
3833 |
|
3834 |
|
3835 |
|
3836 | var sum = jStat.sum;
|
3837 | var range = jStat.arange;
|
3838 |
|
3839 | function qr2(x) {
|
3840 |
|
3841 |
|
3842 |
|
3843 | var n = x.length;
|
3844 | var p = x[0].length;
|
3845 |
|
3846 | var r = jStat.zeros(p, p);
|
3847 | x = jStat.copy(x);
|
3848 |
|
3849 | var i,j,k;
|
3850 | for(j = 0; j < p; j++){
|
3851 | r[j][j] = Math.sqrt(sum(range(n).map(function(i){
|
3852 | return x[i][j] * x[i][j];
|
3853 | })));
|
3854 | for(i = 0; i < n; i++){
|
3855 | x[i][j] = x[i][j] / r[j][j];
|
3856 | }
|
3857 | for(k = j+1; k < p; k++){
|
3858 | r[j][k] = sum(range(n).map(function(i){
|
3859 | return x[i][j] * x[i][k];
|
3860 | }));
|
3861 | for(i = 0; i < n; i++){
|
3862 | x[i][k] = x[i][k] - x[i][j]*r[j][k];
|
3863 | }
|
3864 | }
|
3865 | }
|
3866 | return [x, r];
|
3867 | }
|
3868 |
|
3869 | return qr2;
|
3870 | }()),
|
3871 |
|
3872 | lstsq: (function() {
|
3873 |
|
3874 |
|
3875 |
|
3876 | function R_I(A) {
|
3877 | A = jStat.copy(A);
|
3878 | var size = A.length;
|
3879 | var I = jStat.identity(size);
|
3880 | jStat.arange(size - 1, -1, -1).forEach(function(i) {
|
3881 | jStat.sliceAssign(
|
3882 | I, { row: i }, jStat.divide(jStat.slice(I, { row: i }), A[i][i]));
|
3883 | jStat.sliceAssign(
|
3884 | A, { row: i }, jStat.divide(jStat.slice(A, { row: i }), A[i][i]));
|
3885 | jStat.arange(i).forEach(function(j) {
|
3886 | var c = jStat.multiply(A[j][i], -1);
|
3887 | var Aj = jStat.slice(A, { row: j });
|
3888 | var cAi = jStat.multiply(jStat.slice(A, { row: i }), c);
|
3889 | jStat.sliceAssign(A, { row: j }, jStat.add(Aj, cAi));
|
3890 | var Ij = jStat.slice(I, { row: j });
|
3891 | var cIi = jStat.multiply(jStat.slice(I, { row: i }), c);
|
3892 | jStat.sliceAssign(I, { row: j }, jStat.add(Ij, cIi));
|
3893 | })
|
3894 | });
|
3895 | return I;
|
3896 | }
|
3897 |
|
3898 | function qr_solve(A, b){
|
3899 | var array_mode = false;
|
3900 | if (b[0].length === undefined) {
|
3901 |
|
3902 | b = b.map(function(x){ return [x] });
|
3903 | array_mode = true;
|
3904 | }
|
3905 | var QR = jStat.QR(A);
|
3906 | var Q = QR[0];
|
3907 | var R = QR[1];
|
3908 | var attrs = A[0].length;
|
3909 | var Q1 = jStat.slice(Q,{col:{end:attrs}});
|
3910 | var R1 = jStat.slice(R,{row:{end:attrs}});
|
3911 | var RI = R_I(R1);
|
3912 | var Q2 = jStat.transpose(Q1);
|
3913 |
|
3914 | if(Q2[0].length === undefined){
|
3915 | Q2 = [Q2];
|
3916 | }
|
3917 |
|
3918 | var x = jStat.multiply(jStat.multiply(RI, Q2), b);
|
3919 |
|
3920 | if(x.length === undefined){
|
3921 | x = [[x]];
|
3922 | }
|
3923 |
|
3924 |
|
3925 | if (array_mode)
|
3926 | return x.map(function(i){ return i[0] });
|
3927 | return x;
|
3928 | }
|
3929 |
|
3930 | return qr_solve;
|
3931 | }()),
|
3932 |
|
3933 | jacobi: function jacobi(a) {
|
3934 | var condition = 1;
|
3935 | var n = a.length;
|
3936 | var e = jStat.identity(n, n);
|
3937 | var ev = [];
|
3938 | var b, i, j, p, q, maxim, theta, s;
|
3939 |
|
3940 | while (condition === 1) {
|
3941 | maxim = a[0][1];
|
3942 | p = 0;
|
3943 | q = 1;
|
3944 | for (i = 0; i < n; i++) {
|
3945 | for (j = 0; j < n; j++) {
|
3946 | if (i != j) {
|
3947 | if (maxim < Math.abs(a[i][j])) {
|
3948 | maxim = Math.abs(a[i][j]);
|
3949 | p = i;
|
3950 | q = j;
|
3951 | }
|
3952 | }
|
3953 | }
|
3954 | }
|
3955 | if (a[p][p] === a[q][q])
|
3956 | theta = (a[p][q] > 0) ? Math.PI / 4 : -Math.PI / 4;
|
3957 | else
|
3958 | theta = Math.atan(2 * a[p][q] / (a[p][p] - a[q][q])) / 2;
|
3959 | s = jStat.identity(n, n);
|
3960 | s[p][p] = Math.cos(theta);
|
3961 | s[p][q] = -Math.sin(theta);
|
3962 | s[q][p] = Math.sin(theta);
|
3963 | s[q][q] = Math.cos(theta);
|
3964 |
|
3965 | e = jStat.multiply(e, s);
|
3966 | b = jStat.multiply(jStat.multiply(jStat.inv(s), a), s);
|
3967 | a = b;
|
3968 | condition = 0;
|
3969 | for (i = 1; i < n; i++) {
|
3970 | for (j = 1; j < n; j++) {
|
3971 | if (i != j && Math.abs(a[i][j]) > 0.001) {
|
3972 | condition = 1;
|
3973 | }
|
3974 | }
|
3975 | }
|
3976 | }
|
3977 | for (i = 0; i < n; i++) ev.push(a[i][i]);
|
3978 |
|
3979 | return [e, ev];
|
3980 | },
|
3981 |
|
3982 | rungekutta: function rungekutta(f, h, p, t_j, u_j, order) {
|
3983 | var k1, k2, u_j1, k3, k4;
|
3984 | if (order === 2) {
|
3985 | while (t_j <= p) {
|
3986 | k1 = h * f(t_j, u_j);
|
3987 | k2 = h * f(t_j + h, u_j + k1);
|
3988 | u_j1 = u_j + (k1 + k2) / 2;
|
3989 | u_j = u_j1;
|
3990 | t_j = t_j + h;
|
3991 | }
|
3992 | }
|
3993 | if (order === 4) {
|
3994 | while (t_j <= p) {
|
3995 | k1 = h * f(t_j, u_j);
|
3996 | k2 = h * f(t_j + h / 2, u_j + k1 / 2);
|
3997 | k3 = h * f(t_j + h / 2, u_j + k2 / 2);
|
3998 | k4 = h * f(t_j +h, u_j + k3);
|
3999 | u_j1 = u_j + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
|
4000 | u_j = u_j1;
|
4001 | t_j = t_j + h;
|
4002 | }
|
4003 | }
|
4004 | return u_j;
|
4005 | },
|
4006 |
|
4007 | romberg: function romberg(f, a, b, order) {
|
4008 | var i = 0;
|
4009 | var h = (b - a) / 2;
|
4010 | var x = [];
|
4011 | var h1 = [];
|
4012 | var g = [];
|
4013 | var m, a1, j, k, I;
|
4014 | while (i < order / 2) {
|
4015 | I = f(a);
|
4016 | for (j = a, k = 0; j <= b; j = j + h, k++) x[k] = j;
|
4017 | m = x.length;
|
4018 | for (j = 1; j < m - 1; j++) {
|
4019 | I += (((j % 2) !== 0) ? 4 : 2) * f(x[j]);
|
4020 | }
|
4021 | I = (h / 3) * (I + f(b));
|
4022 | g[i] = I;
|
4023 | h /= 2;
|
4024 | i++;
|
4025 | }
|
4026 | a1 = g.length;
|
4027 | m = 1;
|
4028 | while (a1 !== 1) {
|
4029 | for (j = 0; j < a1 - 1; j++)
|
4030 | h1[j] = ((Math.pow(4, m)) * g[j + 1] - g[j]) / (Math.pow(4, m) - 1);
|
4031 | a1 = h1.length;
|
4032 | g = h1;
|
4033 | h1 = [];
|
4034 | m++;
|
4035 | }
|
4036 | return g;
|
4037 | },
|
4038 |
|
4039 | richardson: function richardson(X, f, x, h) {
|
4040 | function pos(X, x) {
|
4041 | var i = 0;
|
4042 | var n = X.length;
|
4043 | var p;
|
4044 | for (; i < n; i++)
|
4045 | if (X[i] === x) p = i;
|
4046 | return p;
|
4047 | }
|
4048 | var h_min = Math.abs(x - X[pos(X, x) + 1]);
|
4049 | var i = 0;
|
4050 | var g = [];
|
4051 | var h1 = [];
|
4052 | var y1, y2, m, a, j;
|
4053 | while (h >= h_min) {
|
4054 | y1 = pos(X, x + h);
|
4055 | y2 = pos(X, x);
|
4056 | g[i] = (f[y1] - 2 * f[y2] + f[2 * y2 - y1]) / (h * h);
|
4057 | h /= 2;
|
4058 | i++;
|
4059 | }
|
4060 | a = g.length;
|
4061 | m = 1;
|
4062 | while (a != 1) {
|
4063 | for (j = 0; j < a - 1; j++)
|
4064 | h1[j] = ((Math.pow(4, m)) * g[j + 1] - g[j]) / (Math.pow(4, m) - 1);
|
4065 | a = h1.length;
|
4066 | g = h1;
|
4067 | h1 = [];
|
4068 | m++;
|
4069 | }
|
4070 | return g;
|
4071 | },
|
4072 |
|
4073 | simpson: function simpson(f, a, b, n) {
|
4074 | var h = (b - a) / n;
|
4075 | var I = f(a);
|
4076 | var x = [];
|
4077 | var j = a;
|
4078 | var k = 0;
|
4079 | var i = 1;
|
4080 | var m;
|
4081 | for (; j <= b; j = j + h, k++)
|
4082 | x[k] = j;
|
4083 | m = x.length;
|
4084 | for (; i < m - 1; i++) {
|
4085 | I += ((i % 2 !== 0) ? 4 : 2) * f(x[i]);
|
4086 | }
|
4087 | return (h / 3) * (I + f(b));
|
4088 | },
|
4089 |
|
4090 | hermite: function hermite(X, F, dF, value) {
|
4091 | var n = X.length;
|
4092 | var p = 0;
|
4093 | var i = 0;
|
4094 | var l = [];
|
4095 | var dl = [];
|
4096 | var A = [];
|
4097 | var B = [];
|
4098 | var j;
|
4099 | for (; i < n; i++) {
|
4100 | l[i] = 1;
|
4101 | for (j = 0; j < n; j++) {
|
4102 | if (i != j) l[i] *= (value - X[j]) / (X[i] - X[j]);
|
4103 | }
|
4104 | dl[i] = 0;
|
4105 | for (j = 0; j < n; j++) {
|
4106 | if (i != j) dl[i] += 1 / (X [i] - X[j]);
|
4107 | }
|
4108 | A[i] = (1 - 2 * (value - X[i]) * dl[i]) * (l[i] * l[i]);
|
4109 | B[i] = (value - X[i]) * (l[i] * l[i]);
|
4110 | p += (A[i] * F[i] + B[i] * dF[i]);
|
4111 | }
|
4112 | return p;
|
4113 | },
|
4114 |
|
4115 | lagrange: function lagrange(X, F, value) {
|
4116 | var p = 0;
|
4117 | var i = 0;
|
4118 | var j, l;
|
4119 | var n = X.length;
|
4120 | for (; i < n; i++) {
|
4121 | l = F[i];
|
4122 | for (j = 0; j < n; j++) {
|
4123 |
|
4124 | if (i != j) l *= (value - X[j]) / (X[i] - X[j]);
|
4125 | }
|
4126 |
|
4127 | p += l;
|
4128 | }
|
4129 | return p;
|
4130 | },
|
4131 |
|
4132 | cubic_spline: function cubic_spline(X, F, value) {
|
4133 | var n = X.length;
|
4134 | var i = 0, j;
|
4135 | var A = [];
|
4136 | var B = [];
|
4137 | var alpha = [];
|
4138 | var c = [];
|
4139 | var h = [];
|
4140 | var b = [];
|
4141 | var d = [];
|
4142 | for (; i < n - 1; i++)
|
4143 | h[i] = X[i + 1] - X[i];
|
4144 | alpha[0] = 0;
|
4145 | for (i = 1; i < n - 1; i++) {
|
4146 | alpha[i] = (3 / h[i]) * (F[i + 1] - F[i]) -
|
4147 | (3 / h[i-1]) * (F[i] - F[i-1]);
|
4148 | }
|
4149 | for (i = 1; i < n - 1; i++) {
|
4150 | A[i] = [];
|
4151 | B[i] = [];
|
4152 | A[i][i-1] = h[i-1];
|
4153 | A[i][i] = 2 * (h[i - 1] + h[i]);
|
4154 | A[i][i+1] = h[i];
|
4155 | B[i][0] = alpha[i];
|
4156 | }
|
4157 | c = jStat.multiply(jStat.inv(A), B);
|
4158 | for (j = 0; j < n - 1; j++) {
|
4159 | b[j] = (F[j + 1] - F[j]) / h[j] - h[j] * (c[j + 1][0] + 2 * c[j][0]) / 3;
|
4160 | d[j] = (c[j + 1][0] - c[j][0]) / (3 * h[j]);
|
4161 | }
|
4162 | for (j = 0; j < n; j++) {
|
4163 | if (X[j] > value) break;
|
4164 | }
|
4165 | j -= 1;
|
4166 | return F[j] + (value - X[j]) * b[j] + jStat.sq(value-X[j]) *
|
4167 | c[j] + (value - X[j]) * jStat.sq(value - X[j]) * d[j];
|
4168 | },
|
4169 |
|
4170 | gauss_quadrature: function gauss_quadrature() {
|
4171 | throw new Error('gauss_quadrature not yet implemented');
|
4172 | },
|
4173 |
|
4174 | PCA: function PCA(X) {
|
4175 | var m = X.length;
|
4176 | var n = X[0].length;
|
4177 | var i = 0;
|
4178 | var j, temp1;
|
4179 | var u = [];
|
4180 | var D = [];
|
4181 | var result = [];
|
4182 | var temp2 = [];
|
4183 | var Y = [];
|
4184 | var Bt = [];
|
4185 | var B = [];
|
4186 | var C = [];
|
4187 | var V = [];
|
4188 | var Vt = [];
|
4189 | for (i = 0; i < m; i++) {
|
4190 | u[i] = jStat.sum(X[i]) / n;
|
4191 | }
|
4192 | for (i = 0; i < n; i++) {
|
4193 | B[i] = [];
|
4194 | for(j = 0; j < m; j++) {
|
4195 | B[i][j] = X[j][i] - u[j];
|
4196 | }
|
4197 | }
|
4198 | B = jStat.transpose(B);
|
4199 | for (i = 0; i < m; i++) {
|
4200 | C[i] = [];
|
4201 | for (j = 0; j < m; j++) {
|
4202 | C[i][j] = (jStat.dot([B[i]], [B[j]])) / (n - 1);
|
4203 | }
|
4204 | }
|
4205 | result = jStat.jacobi(C);
|
4206 | V = result[0];
|
4207 | D = result[1];
|
4208 | Vt = jStat.transpose(V);
|
4209 | for (i = 0; i < D.length; i++) {
|
4210 | for (j = i; j < D.length; j++) {
|
4211 | if(D[i] < D[j]) {
|
4212 | temp1 = D[i];
|
4213 | D[i] = D[j];
|
4214 | D[j] = temp1;
|
4215 | temp2 = Vt[i];
|
4216 | Vt[i] = Vt[j];
|
4217 | Vt[j] = temp2;
|
4218 | }
|
4219 | }
|
4220 | }
|
4221 | Bt = jStat.transpose(B);
|
4222 | for (i = 0; i < m; i++) {
|
4223 | Y[i] = [];
|
4224 | for (j = 0; j < Bt.length; j++) {
|
4225 | Y[i][j] = jStat.dot([Vt[i]], [Bt[j]]);
|
4226 | }
|
4227 | }
|
4228 | return [X, D, Vt, Y];
|
4229 | }
|
4230 | });
|
4231 |
|
4232 |
|
4233 | (function(funcs) {
|
4234 | for (var i = 0; i < funcs.length; i++) (function(passfunc) {
|
4235 | jStat.fn[passfunc] = function(arg, func) {
|
4236 | var tmpthis = this;
|
4237 |
|
4238 | if (func) {
|
4239 | setTimeout(function() {
|
4240 | func.call(tmpthis, jStat.fn[passfunc].call(tmpthis, arg));
|
4241 | }, 15);
|
4242 | return this;
|
4243 | }
|
4244 | if (typeof jStat[passfunc](this, arg) === 'number')
|
4245 | return jStat[passfunc](this, arg);
|
4246 | else
|
4247 | return jStat(jStat[passfunc](this, arg));
|
4248 | };
|
4249 | }(funcs[i]));
|
4250 | }('add divide multiply subtract dot pow exp log abs norm angle'.split(' ')));
|
4251 |
|
4252 | }(jStat, Math));
|
4253 | (function(jStat, Math) {
|
4254 |
|
4255 | var slice = [].slice;
|
4256 | var isNumber = jStat.utils.isNumber;
|
4257 | var isArray = jStat.utils.isArray;
|
4258 |
|
4259 |
|
4260 |
|
4261 | jStat.extend({
|
4262 |
|
4263 |
|
4264 |
|
4265 | zscore: function zscore() {
|
4266 | var args = slice.call(arguments);
|
4267 | if (isNumber(args[1])) {
|
4268 | return (args[0] - args[1]) / args[2];
|
4269 | }
|
4270 | return (args[0] - jStat.mean(args[1])) / jStat.stdev(args[1], args[2]);
|
4271 | },
|
4272 |
|
4273 |
|
4274 |
|
4275 |
|
4276 |
|
4277 | ztest: function ztest() {
|
4278 | var args = slice.call(arguments);
|
4279 | var z;
|
4280 | if (isArray(args[1])) {
|
4281 |
|
4282 | z = jStat.zscore(args[0],args[1],args[3]);
|
4283 | return (args[2] === 1) ?
|
4284 | (jStat.normal.cdf(-Math.abs(z), 0, 1)) :
|
4285 | (jStat.normal.cdf(-Math.abs(z), 0, 1)*2);
|
4286 | } else {
|
4287 | if (args.length > 2) {
|
4288 |
|
4289 | z = jStat.zscore(args[0],args[1],args[2]);
|
4290 | return (args[3] === 1) ?
|
4291 | (jStat.normal.cdf(-Math.abs(z),0,1)) :
|
4292 | (jStat.normal.cdf(-Math.abs(z),0,1)* 2);
|
4293 | } else {
|
4294 |
|
4295 | z = args[0];
|
4296 | return (args[1] === 1) ?
|
4297 | (jStat.normal.cdf(-Math.abs(z),0,1)) :
|
4298 | (jStat.normal.cdf(-Math.abs(z),0,1)*2);
|
4299 | }
|
4300 | }
|
4301 | }
|
4302 | });
|
4303 |
|
4304 | jStat.extend(jStat.fn, {
|
4305 | zscore: function zscore(value, flag) {
|
4306 | return (value - this.mean()) / this.stdev(flag);
|
4307 | },
|
4308 |
|
4309 | ztest: function ztest(value, sides, flag) {
|
4310 | var zscore = Math.abs(this.zscore(value, flag));
|
4311 | return (sides === 1) ?
|
4312 | (jStat.normal.cdf(-zscore, 0, 1)) :
|
4313 | (jStat.normal.cdf(-zscore, 0, 1) * 2);
|
4314 | }
|
4315 | });
|
4316 |
|
4317 |
|
4318 | jStat.extend({
|
4319 |
|
4320 |
|
4321 |
|
4322 | tscore: function tscore() {
|
4323 | var args = slice.call(arguments);
|
4324 | return (args.length === 4) ?
|
4325 | ((args[0] - args[1]) / (args[2] / Math.sqrt(args[3]))) :
|
4326 | ((args[0] - jStat.mean(args[1])) /
|
4327 | (jStat.stdev(args[1], true) / Math.sqrt(args[1].length)));
|
4328 | },
|
4329 |
|
4330 |
|
4331 |
|
4332 |
|
4333 |
|
4334 | ttest: function ttest() {
|
4335 | var args = slice.call(arguments);
|
4336 | var tscore;
|
4337 | if (args.length === 5) {
|
4338 | tscore = Math.abs(jStat.tscore(args[0], args[1], args[2], args[3]));
|
4339 | return (args[4] === 1) ?
|
4340 | (jStat.studentt.cdf(-tscore, args[3]-1)) :
|
4341 | (jStat.studentt.cdf(-tscore, args[3]-1)*2);
|
4342 | }
|
4343 | if (isNumber(args[1])) {
|
4344 | tscore = Math.abs(args[0])
|
4345 | return (args[2] == 1) ?
|
4346 | (jStat.studentt.cdf(-tscore, args[1]-1)) :
|
4347 | (jStat.studentt.cdf(-tscore, args[1]-1) * 2);
|
4348 | }
|
4349 | tscore = Math.abs(jStat.tscore(args[0], args[1]))
|
4350 | return (args[2] == 1) ?
|
4351 | (jStat.studentt.cdf(-tscore, args[1].length-1)) :
|
4352 | (jStat.studentt.cdf(-tscore, args[1].length-1) * 2);
|
4353 | }
|
4354 | });
|
4355 |
|
4356 | jStat.extend(jStat.fn, {
|
4357 | tscore: function tscore(value) {
|
4358 | return (value - this.mean()) / (this.stdev(true) / Math.sqrt(this.cols()));
|
4359 | },
|
4360 |
|
4361 | ttest: function ttest(value, sides) {
|
4362 | return (sides === 1) ?
|
4363 | (1 - jStat.studentt.cdf(Math.abs(this.tscore(value)), this.cols()-1)) :
|
4364 | (jStat.studentt.cdf(-Math.abs(this.tscore(value)), this.cols()-1)*2);
|
4365 | }
|
4366 | });
|
4367 |
|
4368 |
|
4369 | jStat.extend({
|
4370 |
|
4371 |
|
4372 |
|
4373 |
|
4374 | anovafscore: function anovafscore() {
|
4375 | var args = slice.call(arguments),
|
4376 | expVar, sample, sampMean, sampSampMean, tmpargs, unexpVar, i, j;
|
4377 | if (args.length === 1) {
|
4378 | tmpargs = new Array(args[0].length);
|
4379 | for (i = 0; i < args[0].length; i++) {
|
4380 | tmpargs[i] = args[0][i];
|
4381 | }
|
4382 | args = tmpargs;
|
4383 | }
|
4384 |
|
4385 | sample = new Array();
|
4386 | for (i = 0; i < args.length; i++) {
|
4387 | sample = sample.concat(args[i]);
|
4388 | }
|
4389 | sampMean = jStat.mean(sample);
|
4390 |
|
4391 | expVar = 0;
|
4392 | for (i = 0; i < args.length; i++) {
|
4393 | expVar = expVar + args[i].length * Math.pow(jStat.mean(args[i]) - sampMean, 2);
|
4394 | }
|
4395 | expVar /= (args.length - 1);
|
4396 |
|
4397 | unexpVar = 0;
|
4398 | for (i = 0; i < args.length; i++) {
|
4399 | sampSampMean = jStat.mean(args[i]);
|
4400 | for (j = 0; j < args[i].length; j++) {
|
4401 | unexpVar += Math.pow(args[i][j] - sampSampMean, 2);
|
4402 | }
|
4403 | }
|
4404 | unexpVar /= (sample.length - args.length);
|
4405 | return expVar / unexpVar;
|
4406 | },
|
4407 |
|
4408 |
|
4409 |
|
4410 |
|
4411 | anovaftest: function anovaftest() {
|
4412 | var args = slice.call(arguments),
|
4413 | df1, df2, n, i;
|
4414 | if (isNumber(args[0])) {
|
4415 | return 1 - jStat.centralF.cdf(args[0], args[1], args[2]);
|
4416 | }
|
4417 | var anovafscore = jStat.anovafscore(args);
|
4418 | df1 = args.length - 1;
|
4419 | n = 0;
|
4420 | for (i = 0; i < args.length; i++) {
|
4421 | n = n + args[i].length;
|
4422 | }
|
4423 | df2 = n - df1 - 1;
|
4424 | return 1 - jStat.centralF.cdf(anovafscore, df1, df2);
|
4425 | },
|
4426 |
|
4427 | ftest: function ftest(fscore, df1, df2) {
|
4428 | return 1 - jStat.centralF.cdf(fscore, df1, df2);
|
4429 | }
|
4430 | });
|
4431 |
|
4432 | jStat.extend(jStat.fn, {
|
4433 | anovafscore: function anovafscore() {
|
4434 | return jStat.anovafscore(this.toArray());
|
4435 | },
|
4436 |
|
4437 | anovaftes: function anovaftes() {
|
4438 | var n = 0;
|
4439 | var i;
|
4440 | for (i = 0; i < this.length; i++) {
|
4441 | n = n + this[i].length;
|
4442 | }
|
4443 | return jStat.ftest(this.anovafscore(), this.length - 1, n - this.length);
|
4444 | }
|
4445 | });
|
4446 |
|
4447 |
|
4448 | jStat.extend({
|
4449 |
|
4450 |
|
4451 |
|
4452 | qscore: function qscore() {
|
4453 | var args = slice.call(arguments);
|
4454 | var mean1, mean2, n1, n2, sd;
|
4455 | if (isNumber(args[0])) {
|
4456 | mean1 = args[0];
|
4457 | mean2 = args[1];
|
4458 | n1 = args[2];
|
4459 | n2 = args[3];
|
4460 | sd = args[4];
|
4461 | } else {
|
4462 | mean1 = jStat.mean(args[0]);
|
4463 | mean2 = jStat.mean(args[1]);
|
4464 | n1 = args[0].length;
|
4465 | n2 = args[1].length;
|
4466 | sd = args[2];
|
4467 | }
|
4468 | return Math.abs(mean1 - mean2) / (sd * Math.sqrt((1 / n1 + 1 / n2) / 2));
|
4469 | },
|
4470 |
|
4471 |
|
4472 |
|
4473 |
|
4474 |
|
4475 | qtest: function qtest() {
|
4476 | var args = slice.call(arguments);
|
4477 |
|
4478 | var qscore;
|
4479 | if (args.length === 3) {
|
4480 | qscore = args[0];
|
4481 | args = args.slice(1);
|
4482 | } else if (args.length === 7) {
|
4483 | qscore = jStat.qscore(args[0], args[1], args[2], args[3], args[4]);
|
4484 | args = args.slice(5);
|
4485 | } else {
|
4486 | qscore = jStat.qscore(args[0], args[1], args[2]);
|
4487 | args = args.slice(3);
|
4488 | }
|
4489 |
|
4490 | var n = args[0];
|
4491 | var k = args[1];
|
4492 |
|
4493 | return 1 - jStat.tukey.cdf(qscore, k, n - k);
|
4494 | },
|
4495 |
|
4496 | tukeyhsd: function tukeyhsd(arrays) {
|
4497 | var sd = jStat.pooledstdev(arrays);
|
4498 | var means = arrays.map(function (arr) {return jStat.mean(arr);});
|
4499 | var n = arrays.reduce(function (n, arr) {return n + arr.length;}, 0);
|
4500 |
|
4501 | var results = [];
|
4502 | for (var i = 0; i < arrays.length; ++i) {
|
4503 | for (var j = i + 1; j < arrays.length; ++j) {
|
4504 | var p = jStat.qtest(means[i], means[j], arrays[i].length, arrays[j].length, sd, n, arrays.length);
|
4505 | results.push([[i, j], p]);
|
4506 | }
|
4507 | }
|
4508 |
|
4509 | return results;
|
4510 | }
|
4511 | });
|
4512 |
|
4513 |
|
4514 | jStat.extend({
|
4515 |
|
4516 |
|
4517 |
|
4518 | normalci: function normalci() {
|
4519 | var args = slice.call(arguments),
|
4520 | ans = new Array(2),
|
4521 | change;
|
4522 | if (args.length === 4) {
|
4523 | change = Math.abs(jStat.normal.inv(args[1] / 2, 0, 1) *
|
4524 | args[2] / Math.sqrt(args[3]));
|
4525 | } else {
|
4526 | change = Math.abs(jStat.normal.inv(args[1] / 2, 0, 1) *
|
4527 | jStat.stdev(args[2]) / Math.sqrt(args[2].length));
|
4528 | }
|
4529 | ans[0] = args[0] - change;
|
4530 | ans[1] = args[0] + change;
|
4531 | return ans;
|
4532 | },
|
4533 |
|
4534 |
|
4535 |
|
4536 |
|
4537 | tci: function tci() {
|
4538 | var args = slice.call(arguments),
|
4539 | ans = new Array(2),
|
4540 | change;
|
4541 | if (args.length === 4) {
|
4542 | change = Math.abs(jStat.studentt.inv(args[1] / 2, args[3] - 1) *
|
4543 | args[2] / Math.sqrt(args[3]));
|
4544 | } else {
|
4545 | change = Math.abs(jStat.studentt.inv(args[1] / 2, args[2].length - 1) *
|
4546 | jStat.stdev(args[2], true) / Math.sqrt(args[2].length));
|
4547 | }
|
4548 | ans[0] = args[0] - change;
|
4549 | ans[1] = args[0] + change;
|
4550 | return ans;
|
4551 | },
|
4552 |
|
4553 | significant: function significant(pvalue, alpha) {
|
4554 | return pvalue < alpha;
|
4555 | }
|
4556 | });
|
4557 |
|
4558 | jStat.extend(jStat.fn, {
|
4559 | normalci: function normalci(value, alpha) {
|
4560 | return jStat.normalci(value, alpha, this.toArray());
|
4561 | },
|
4562 |
|
4563 | tci: function tci(value, alpha) {
|
4564 | return jStat.tci(value, alpha, this.toArray());
|
4565 | }
|
4566 | });
|
4567 |
|
4568 |
|
4569 | function differenceOfProportions(p1, n1, p2, n2) {
|
4570 | if (p1 > 1 || p2 > 1 || p1 <= 0 || p2 <= 0) {
|
4571 | throw new Error("Proportions should be greater than 0 and less than 1")
|
4572 | }
|
4573 | var pooled = (p1 * n1 + p2 * n2) / (n1 + n2);
|
4574 | var se = Math.sqrt(pooled * (1 - pooled) * ((1/n1) + (1/n2)));
|
4575 | return (p1 - p2) / se;
|
4576 | }
|
4577 |
|
4578 |
|
4579 | jStat.extend(jStat.fn, {
|
4580 | oneSidedDifferenceOfProportions: function oneSidedDifferenceOfProportions(p1, n1, p2, n2) {
|
4581 | var z = differenceOfProportions(p1, n1, p2, n2);
|
4582 | return jStat.ztest(z, 1);
|
4583 | },
|
4584 |
|
4585 | twoSidedDifferenceOfProportions: function twoSidedDifferenceOfProportions(p1, n1, p2, n2) {
|
4586 | var z = differenceOfProportions(p1, n1, p2, n2);
|
4587 | return jStat.ztest(z, 2);
|
4588 | }
|
4589 | });
|
4590 |
|
4591 | }(jStat, Math));
|
4592 | jStat.models = (function(){
|
4593 | function sub_regress(exog) {
|
4594 | var var_count = exog[0].length;
|
4595 | var modelList = jStat.arange(var_count).map(function(endog_index) {
|
4596 | var exog_index =
|
4597 | jStat.arange(var_count).filter(function(i){return i!==endog_index});
|
4598 | return ols(jStat.col(exog, endog_index).map(function(x){ return x[0] }),
|
4599 | jStat.col(exog, exog_index))
|
4600 | });
|
4601 | return modelList;
|
4602 | }
|
4603 |
|
4604 |
|
4605 |
|
4606 |
|
4607 |
|
4608 | function ols(endog, exog) {
|
4609 | var nobs = endog.length;
|
4610 | var df_model = exog[0].length - 1;
|
4611 | var df_resid = nobs-df_model - 1;
|
4612 | var coef = jStat.lstsq(exog, endog);
|
4613 | var predict =
|
4614 | jStat.multiply(exog, coef.map(function(x) { return [x] }))
|
4615 | .map(function(p) { return p[0] });
|
4616 | var resid = jStat.subtract(endog, predict);
|
4617 | var ybar = jStat.mean(endog);
|
4618 |
|
4619 |
|
4620 |
|
4621 |
|
4622 | var SSE = jStat.sum(predict.map(function(f) {
|
4623 | return Math.pow(f - ybar, 2);
|
4624 | }));
|
4625 | var SSR = jStat.sum(endog.map(function(y, i) {
|
4626 | return Math.pow(y - predict[i], 2);
|
4627 | }));
|
4628 | var SST = SSE + SSR;
|
4629 | var R2 = (SSE / SST);
|
4630 | return {
|
4631 | exog:exog,
|
4632 | endog:endog,
|
4633 | nobs:nobs,
|
4634 | df_model:df_model,
|
4635 | df_resid:df_resid,
|
4636 | coef:coef,
|
4637 | predict:predict,
|
4638 | resid:resid,
|
4639 | ybar:ybar,
|
4640 | SST:SST,
|
4641 | SSE:SSE,
|
4642 | SSR:SSR,
|
4643 | R2:R2
|
4644 | };
|
4645 | }
|
4646 |
|
4647 |
|
4648 |
|
4649 | function t_test(model) {
|
4650 | var subModelList = sub_regress(model.exog);
|
4651 |
|
4652 | var sigmaHat = Math.sqrt(model.SSR / (model.df_resid));
|
4653 | var seBetaHat = subModelList.map(function(mod) {
|
4654 | var SST = mod.SST;
|
4655 | var R2 = mod.R2;
|
4656 | return sigmaHat / Math.sqrt(SST * (1 - R2));
|
4657 | });
|
4658 | var tStatistic = model.coef.map(function(coef, i) {
|
4659 | return (coef - 0) / seBetaHat[i];
|
4660 | });
|
4661 | var pValue = tStatistic.map(function(t) {
|
4662 | var leftppf = jStat.studentt.cdf(t, model.df_resid);
|
4663 | return (leftppf > 0.5 ? 1 - leftppf : leftppf) * 2;
|
4664 | });
|
4665 | var c = jStat.studentt.inv(0.975, model.df_resid);
|
4666 | var interval95 = model.coef.map(function(coef, i) {
|
4667 | var d = c * seBetaHat[i];
|
4668 | return [coef - d, coef + d];
|
4669 | })
|
4670 | return {
|
4671 | se: seBetaHat,
|
4672 | t: tStatistic,
|
4673 | p: pValue,
|
4674 | sigmaHat: sigmaHat,
|
4675 | interval95: interval95
|
4676 | };
|
4677 | }
|
4678 |
|
4679 | function F_test(model) {
|
4680 | var F_statistic =
|
4681 | (model.R2 / model.df_model) / ((1 - model.R2) / model.df_resid);
|
4682 | var fcdf = function(x, n1, n2) {
|
4683 | return jStat.beta.cdf(x / (n2 / n1 + x), n1 / 2, n2 / 2)
|
4684 | }
|
4685 | var pvalue = 1 - fcdf(F_statistic, model.df_model, model.df_resid);
|
4686 | return { F_statistic: F_statistic, pvalue: pvalue };
|
4687 | }
|
4688 |
|
4689 | function ols_wrap(endog, exog) {
|
4690 | var model = ols(endog,exog);
|
4691 | var ttest = t_test(model);
|
4692 | var ftest = F_test(model);
|
4693 |
|
4694 |
|
4695 | var adjust_R2 =
|
4696 | 1 - (1 - model.R2) * ((model.nobs - 1) / (model.df_resid));
|
4697 | model.t = ttest;
|
4698 | model.f = ftest;
|
4699 | model.adjust_R2 = adjust_R2;
|
4700 | return model;
|
4701 | }
|
4702 |
|
4703 | return { ols: ols_wrap };
|
4704 | })();
|
4705 |
|
4706 |
|
4707 |
|
4708 |
|
4709 |
|
4710 |
|
4711 |
|
4712 |
|
4713 |
|
4714 |
|
4715 |
|
4716 | jStat.extend({
|
4717 | buildxmatrix: function buildxmatrix(){
|
4718 |
|
4719 |
|
4720 |
|
4721 |
|
4722 | var matrixRows = new Array(arguments.length);
|
4723 | for(var i=0;i<arguments.length;i++){
|
4724 | var array = [1];
|
4725 | matrixRows[i]= array.concat(arguments[i]);
|
4726 | }
|
4727 | return jStat(matrixRows);
|
4728 |
|
4729 | },
|
4730 |
|
4731 | builddxmatrix: function builddxmatrix() {
|
4732 |
|
4733 |
|
4734 | var matrixRows = new Array(arguments[0].length);
|
4735 | for(var i=0;i<arguments[0].length;i++){
|
4736 | var array = [1]
|
4737 | matrixRows[i]= array.concat(arguments[0][i]);
|
4738 | }
|
4739 | return jStat(matrixRows);
|
4740 |
|
4741 | },
|
4742 |
|
4743 | buildjxmatrix: function buildjxmatrix(jMat) {
|
4744 |
|
4745 | var pass = new Array(jMat.length)
|
4746 | for(var i=0;i<jMat.length;i++){
|
4747 | pass[i] = jMat[i];
|
4748 | }
|
4749 | return jStat.builddxmatrix(pass);
|
4750 |
|
4751 | },
|
4752 |
|
4753 | buildymatrix: function buildymatrix(array){
|
4754 | return jStat(array).transpose();
|
4755 | },
|
4756 |
|
4757 | buildjymatrix: function buildjymatrix(jMat){
|
4758 | return jMat.transpose();
|
4759 | },
|
4760 |
|
4761 | matrixmult: function matrixmult(A,B){
|
4762 | var i, j, k, result, sum;
|
4763 | if (A.cols() == B.rows()) {
|
4764 | if(B.rows()>1){
|
4765 | result = [];
|
4766 | for (i = 0; i < A.rows(); i++) {
|
4767 | result[i] = [];
|
4768 | for (j = 0; j < B.cols(); j++) {
|
4769 | sum = 0;
|
4770 | for (k = 0; k < A.cols(); k++) {
|
4771 | sum += A.toArray()[i][k] * B.toArray()[k][j];
|
4772 | }
|
4773 | result[i][j] = sum;
|
4774 | }
|
4775 | }
|
4776 | return jStat(result);
|
4777 | }
|
4778 | result = [];
|
4779 | for (i = 0; i < A.rows(); i++) {
|
4780 | result[i] = [];
|
4781 | for (j = 0; j < B.cols(); j++) {
|
4782 | sum = 0;
|
4783 | for (k = 0; k < A.cols(); k++) {
|
4784 | sum += A.toArray()[i][k] * B.toArray()[j];
|
4785 | }
|
4786 | result[i][j] = sum;
|
4787 | }
|
4788 | }
|
4789 | return jStat(result);
|
4790 | }
|
4791 | },
|
4792 |
|
4793 |
|
4794 |
|
4795 | regress: function regress(jMatX,jMatY){
|
4796 |
|
4797 |
|
4798 | var innerinv = jStat.xtranspxinv(jMatX);
|
4799 |
|
4800 | var xtransp = jMatX.transpose();
|
4801 | var next = jStat.matrixmult(jStat(innerinv),xtransp);
|
4802 | return jStat.matrixmult(next,jMatY);
|
4803 |
|
4804 | },
|
4805 |
|
4806 | regresst: function regresst(jMatX,jMatY,sides){
|
4807 | var beta = jStat.regress(jMatX,jMatY);
|
4808 |
|
4809 | var compile = {};
|
4810 | compile.anova = {};
|
4811 | var jMatYBar = jStat.jMatYBar(jMatX, beta);
|
4812 | compile.yBar = jMatYBar;
|
4813 | var yAverage = jMatY.mean();
|
4814 | compile.anova.residuals = jStat.residuals(jMatY, jMatYBar);
|
4815 |
|
4816 | compile.anova.ssr = jStat.ssr(jMatYBar, yAverage);
|
4817 | compile.anova.msr = compile.anova.ssr / (jMatX[0].length - 1);
|
4818 |
|
4819 | compile.anova.sse = jStat.sse(jMatY, jMatYBar);
|
4820 | compile.anova.mse =
|
4821 | compile.anova.sse / (jMatY.length - (jMatX[0].length - 1) - 1);
|
4822 |
|
4823 | compile.anova.sst = jStat.sst(jMatY, yAverage);
|
4824 | compile.anova.mst = compile.anova.sst / (jMatY.length - 1);
|
4825 |
|
4826 | compile.anova.r2 = 1 - (compile.anova.sse / compile.anova.sst);
|
4827 | if (compile.anova.r2 < 0) compile.anova.r2 = 0;
|
4828 |
|
4829 | compile.anova.fratio = compile.anova.msr / compile.anova.mse;
|
4830 | compile.anova.pvalue =
|
4831 | jStat.anovaftest(compile.anova.fratio,
|
4832 | jMatX[0].length - 1,
|
4833 | jMatY.length - (jMatX[0].length - 1) - 1);
|
4834 |
|
4835 | compile.anova.rmse = Math.sqrt(compile.anova.mse);
|
4836 |
|
4837 | compile.anova.r2adj = 1 - (compile.anova.mse / compile.anova.mst);
|
4838 | if (compile.anova.r2adj < 0) compile.anova.r2adj = 0;
|
4839 |
|
4840 | compile.stats = new Array(jMatX[0].length);
|
4841 | var covar = jStat.xtranspxinv(jMatX);
|
4842 | var sds, ts, ps;
|
4843 |
|
4844 | for(var i=0; i<beta.length;i++){
|
4845 | sds=Math.sqrt(compile.anova.mse * Math.abs(covar[i][i]));
|
4846 | ts= Math.abs(beta[i] / sds);
|
4847 | ps= jStat.ttest(ts, jMatY.length - jMatX[0].length - 1, sides);
|
4848 |
|
4849 | compile.stats[i]=[beta[i], sds, ts, ps];
|
4850 | }
|
4851 |
|
4852 | compile.regress = beta;
|
4853 | return compile;
|
4854 | },
|
4855 |
|
4856 | xtranspx: function xtranspx(jMatX){
|
4857 | return jStat.matrixmult(jMatX.transpose(),jMatX);
|
4858 | },
|
4859 |
|
4860 |
|
4861 | xtranspxinv: function xtranspxinv(jMatX){
|
4862 | var inner = jStat.matrixmult(jMatX.transpose(),jMatX);
|
4863 | var innerinv = jStat.inv(inner);
|
4864 | return innerinv;
|
4865 | },
|
4866 |
|
4867 | jMatYBar: function jMatYBar(jMatX, beta) {
|
4868 | var yBar = jStat.matrixmult(jMatX, beta);
|
4869 | return new jStat(yBar);
|
4870 | },
|
4871 |
|
4872 | residuals: function residuals(jMatY, jMatYBar) {
|
4873 | return jStat.matrixsubtract(jMatY, jMatYBar);
|
4874 | },
|
4875 |
|
4876 | ssr: function ssr(jMatYBar, yAverage) {
|
4877 | var ssr = 0;
|
4878 | for(var i = 0; i < jMatYBar.length; i++) {
|
4879 | ssr += Math.pow(jMatYBar[i] - yAverage, 2);
|
4880 | }
|
4881 | return ssr;
|
4882 | },
|
4883 |
|
4884 | sse: function sse(jMatY, jMatYBar) {
|
4885 | var sse = 0;
|
4886 | for(var i = 0; i < jMatY.length; i++) {
|
4887 | sse += Math.pow(jMatY[i] - jMatYBar[i], 2);
|
4888 | }
|
4889 | return sse;
|
4890 | },
|
4891 |
|
4892 | sst: function sst(jMatY, yAverage) {
|
4893 | var sst = 0;
|
4894 | for(var i = 0; i < jMatY.length; i++) {
|
4895 | sst += Math.pow(jMatY[i] - yAverage, 2);
|
4896 | }
|
4897 | return sst;
|
4898 | },
|
4899 |
|
4900 | matrixsubtract: function matrixsubtract(A,B){
|
4901 | var ans = new Array(A.length);
|
4902 | for(var i=0;i<A.length;i++){
|
4903 | ans[i] = new Array(A[i].length);
|
4904 | for(var j=0;j<A[i].length;j++){
|
4905 | ans[i][j]=A[i][j]-B[i][j];
|
4906 | }
|
4907 | }
|
4908 | return jStat(ans);
|
4909 | }
|
4910 | });
|
4911 |
|
4912 | jStat.jStat = jStat;
|
4913 |
|
4914 | return jStat;
|
4915 | });
|