Mercurial > repos > adam-novak > hexagram

comparison hexagram-6ae12361157c/hexagram/jstat-1.0.0.js @ 0:1407e3634bcf draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded r11 from test tool shed.
author adam-novak
date Tue, 22 Oct 2013 14:17:59 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:1407e3634bcf
1 function jstat(){}
2 j = jstat;
3 /* Simple JavaScript Inheritance
4 * By John Resig http://ejohn.org/
5 * MIT Licensed.
6 */
7 // Inspired by base2 and Prototype
8 (function(){
9 var initializing = false, fnTest = /xyz/.test(function(){
10 xyz;
11 }) ? /\b_super\b/ : /.*/;
12 // The base Class implementation (does nothing)
13 this.Class = function(){};
14
15 // Create a new Class that inherits from this class
16 Class.extend = function(prop) {
17 var _super = this.prototype;
18
19 // Instantiate a base class (but only create the instance,
20 // don't run the init constructor)
21 initializing = true;
22 var prototype = new this();
23 initializing = false;
24
25 // Copy the properties over onto the new prototype
26 for (var name in prop) {
27 // Check if we're overwriting an existing function
28 prototype[name] = typeof prop[name] == "function" &&
29 typeof _super[name] == "function" && fnTest.test(prop[name]) ?
30 (function(name, fn){
31 return function() {
32 var tmp = this._super;
33
34 // Add a new ._super() method that is the same method
35 // but on the super-class
36 this._super = _super[name];
37
38 // The method only need to be bound temporarily, so we
39 // remove it when we're done executing
40 var ret = fn.apply(this, arguments);
41 this._super = tmp;
42
43 return ret;
44 };
45 })(name, prop[name]) :
46 prop[name];
47 }
48
49 // The dummy class constructor
50 function Class() {
51 // All construction is actually done in the init method
52 if ( !initializing && this.init )
53 this.init.apply(this, arguments);
54 }
55
56 // Populate our constructed prototype object
57 Class.prototype = prototype;
58
59 // Enforce the constructor to be what we expect
60 Class.constructor = Class;
61
62 // And make this class extendable
63 Class.extend = arguments.callee;
64
65 return Class;
66 };
67 })();
68
69 /******************************************************************************/
70 /* Constants */
71 /******************************************************************************/
72 jstat.ONE_SQRT_2PI = 0.3989422804014327;
73 jstat.LN_SQRT_2PI = 0.9189385332046727417803297;
74 jstat.LN_SQRT_PId2 = 0.225791352644727432363097614947;
75 jstat.DBL_MIN = 2.22507e-308;
76 jstat.DBL_EPSILON = 2.220446049250313e-16;
77 jstat.SQRT_32 = 5.656854249492380195206754896838;
78 jstat.TWO_PI = 6.283185307179586;
79 jstat.DBL_MIN_EXP = -999;
80 jstat.SQRT_2dPI = 0.79788456080287;
81 jstat.LN_SQRT_PI = 0.5723649429247;
82 /******************************************************************************/
83 /* jstat Functions */
84 /******************************************************************************/
85 jstat.seq = function(min, max, length) {
86 var r = new Range(min, max, length);
87 return r.getPoints();
88 }
89
90 jstat.dnorm = function(x, mean, sd, log) {
91 if(mean == null) mean = 0;
92 if(sd == null) sd = 1;
93 if(log == null) log = false;
94 var n = new NormalDistribution(mean, sd);
95 if(!isNaN(x)) {
96 // is a number
97 return n._pdf(x, log);
98 } else if(x.length) {
99 var res = [];
100 for(var i = 0; i < x.length; i++) {
101 res.push(n._pdf(x[i], log));
102 }
103 return res;
104 } else {
105 throw "Illegal argument: x";
106 }
107 }
108
109 jstat.pnorm = function(q, mean, sd, lower_tail, log) {
110 if(mean == null) mean = 0;
111 if(sd == null) sd = 1;
112 if(lower_tail == null) lower_tail = true;
113 if(log == null) log = false;
114
115 var n = new NormalDistribution(mean, sd);
116 if(!isNaN(q)) {
117 // is a number
118 return n._cdf(q, lower_tail, log);
119 } else if(q.length) {
120 var res = [];
121 for(var i = 0; i < q.length; i++) {
122 res.push(n._cdf(q[i], lower_tail, log));
123 }
124 return res;
125 } else {
126 throw "Illegal argument: x";
127 }
128 }
129
130 jstat.dlnorm = function(x, meanlog, sdlog, log) {
131 if(meanlog == null) meanlog = 0;
132 if(sdlog == null) sdlog = 1;
133 if(log == null) log = false;
134 var n = new LogNormalDistribution(meanlog, sdlog);
135 if(!isNaN(x)) {
136 // is a number
137 return n._pdf(x, log);
138 } else if(x.length) {
139 var res = [];
140 for(var i = 0; i < x.length; i++) {
141 res.push(n._pdf(x[i], log));
142 }
143 return res;
144 } else {
145 throw "Illegal argument: x";
146 }
147 }
148
149 jstat.plnorm = function(q, meanlog, sdlog, lower_tail, log) {
150 if(meanlog == null) meanlog = 0;
151 if(sdlog == null) sdlog = 1;
152 if(lower_tail == null) lower_tail = true;
153 if(log == null) log = false;
154
155 var n = new LogNormalDistribution(meanlog, sdlog);
156 if(!isNaN(q)) {
157 // is a number
158 return n._cdf(q, lower_tail, log);
159 }
160 else if(q.length) {
161 var res = [];
162 for(var i = 0; i < q.length; i++) {
163 res.push(n._cdf(q[i], lower_tail, log));
164 }
165 return res;
166 } else {
167 throw "Illegal argument: x";
168 }
169 }
170
171 jstat.dbeta = function(x, alpha, beta, ncp, log) {
172 if(ncp == null) ncp = 0;
173 if(log == null) log = false;
174 var b = new BetaDistribution(alpha, beta);
175 if(!isNaN(x)) {
176 // is a number
177 return b._pdf(x, log);
178 }
179 else if(x.length) {
180 var res = [];
181 for(var i = 0; i < x.length; i++) {
182 res.push(b._pdf(x[i], log));
183 }
184 return res;
185 } else {
186 throw "Illegal argument: x";
187 }
188 }
189
190 jstat.pbeta = function(q, alpha, beta, ncp, lower_tail, log) {
191 if(ncp == null) ncp = 0;
192 if(log == null) log = false;
193 if(lower_tail == null) lower_tail = true;
194
195 var b = new BetaDistribution(alpha, beta);
196 if(!isNaN(q)) {
197 // is a number
198 return b._cdf(q, lower_tail, log);
199 } else if(q.length) {
200 var res = [];
201 for(var i = 0; i < q.length; i++) {
202 res.push(b._cdf(q[i], lower_tail, log));
203 }
204 return res;
205 }
206 else {
207 throw "Illegal argument: x";
208 }
209 }
210
211 jstat.dgamma = function(x, shape, rate, scale, log) {
212 if(rate == null) rate = 1;
213 if(scale == null) scale = 1/rate;
214 if(log == null) log = false;
215
216 var g = new GammaDistribution(shape, scale);
217 if(!isNaN(x)) {
218 // is a number
219 return g._pdf(x, log);
220 } else if(x.length) {
221 var res = [];
222 for(var i = 0; i < x.length; i++) {
223 res.push(g._pdf(x[i], log));
224 }
225 return res;
226 } else {
227 throw "Illegal argument: x";
228 }
229 }
230
231 jstat.pgamma = function(q, shape, rate, scale, lower_tail, log) {
232 if(rate == null) rate = 1;
233 if(scale == null) scale = 1/rate;
234 if(lower_tail == null) lower_tail = true;
235 if(log == null) log = false;
236
237 var g = new GammaDistribution(shape, scale);
238 if(!isNaN(q)) {
239 // is a number
240 return g._cdf(q, lower_tail, log);
241 } else if(q.length) {
242 var res = [];
243 for(var i = 0; i < q.length; i++) {
244 res.push(g._cdf(q[i], lower_tail, log));
245 }
246 return res;
247 } else {
248 throw "Illegal argument: x";
249 }
250
251 }
252
253 jstat.dt = function(x, df, ncp, log) {
254 if(log == null) log = false;
255
256 var t = new StudentTDistribution(df, ncp);
257 if(!isNaN(x)) {
258 // is a number
259 return t._pdf(x, log);
260 } else if(x.length) {
261 var res = [];
262 for(var i = 0; i < x.length; i++) {
263 res.push(t._pdf(x[i], log));
264 }
265 return res;
266 } else {
267 throw "Illegal argument: x";
268 }
269
270 }
271
272 jstat.pt = function(q, df, ncp, lower_tail, log) {
273 if(lower_tail == null) lower_tail = true;
274 if(log == null) log = false;
275
276 var t = new StudentTDistribution(df, ncp);
277 if(!isNaN(q)) {
278 // is a number
279 return t._cdf(q, lower_tail, log);
280 } else if(q.length) {
281 var res = [];
282 for(var i = 0; i < q.length; i++) {
283 res.push(t._cdf(q[i], lower_tail, log));
284 }
285 return res;
286 } else {
287 throw "Illegal argument: x";
288 }
289
290 }
291
292 jstat.plot = function(x, y, options) {
293 if(x == null) {
294 throw "x is undefined in jstat.plot";
295 }
296 if(y == null) {
297 throw "y is undefined in jstat.plot";
298 }
299 if(x.length != y.length) {
300 throw "x and y lengths differ in jstat.plot";
301 }
302
303 var flotOpt = {
304 series: {
305 lines: {
306
307 },
308 points: {
309
310 }
311 }
312 };
313
314 // combine x & y
315 var series = [];
316 if(x.length == undefined) {
317 // single point
318 series.push([x, y]);
319 flotOpt.series.points.show = true;
320 } else {
321 // array
322 for(var i = 0; i < x.length; i++) {
323 series.push([x[i], y[i]]);
324 }
325 }
326
327 var title = 'jstat graph';
328
329 // configure Flot options
330 if(options != null) {
331 // options = JSON.parse(String(options));
332 if(options.type != null) {
333 if(options.type == 'l') {
334 flotOpt.series.lines.show = true;
335 } else if (options.type == 'p') {
336 flotOpt.series.lines.show = false;
337 flotOpt.series.points.show = true;
338 }
339 }
340 if(options.hover != null) {
341 flotOpt.grid = {
342 hoverable: options.hover
343 }
344 }
345
346 if(options.main != null) {
347 title = options.main;
348 }
349 }
350 var now = new Date();
351 var hash = now.getMilliseconds() * now.getMinutes() + now.getSeconds();
352 $('body').append('<div title="' + title + '" style="display: none;" id="'+ hash +'"><div id="graph-' + hash + '" style="width:95%; height: 95%"></div></div>');
353
354 $('#' + hash).dialog({
355 modal: false,
356 width: 475,
357 height: 475,
358 resizable: true,
359 resize: function() {
360 $.plot($('#graph-' + hash), [series], flotOpt);
361 },
362 open: function(event, ui) {
363 var id = '#graph-' + hash;
364 $.plot($('#graph-' + hash), [series], flotOpt);
365 }
366 })
367 }
368
369 /******************************************************************************/
370 /* Special Functions */
371 /******************************************************************************/
372
373 jstat.log10 = function(arg) {
374 return Math.log(arg) / Math.LN10;
375 }
376
377 /*
378 *
379 */
380 jstat.toSigFig = function(num, n) {
381 if(num == 0) {
382 return 0;
383 }
384 var d = Math.ceil(jstat.log10(num < 0 ? -num: num));
385 var power = n - parseInt(d);
386 var magnitude = Math.pow(10,power);
387 var shifted = Math.round(num*magnitude);
388 return shifted/magnitude;
389 }
390
391 jstat.trunc = function(x) {
392 return (x > 0) ? Math.floor(x) : Math.ceil(x);
393 }
394
395 /**
396 * Tests whether x is a finite number
397 */
398 jstat.isFinite = function(x) {
399 return (!isNaN(x) && (x != Number.POSITIVE_INFINITY) && (x != Number.NEGATIVE_INFINITY));
400 }
401
402 /**
403 * dopois_raw() computes the Poisson probability lb^x exp(-lb) / x!.
404 * This does not check that x is an integer, since dgamma() may
405 * call this with a fractional x argument. Any necessary argument
406 * checks should be done in the calling function.
407 */
408 jstat.dopois_raw = function(x, lambda, give_log) {
409 /* x >= 0 ; integer for dpois(), but not e.g. for pgamma()!
410 lambda >= 0
411 */
412 if (lambda == 0) {
413 if(x == 0) {
414 return(give_log) ? 0.0 : 1.0; //R_D__1
415 }
416 return (give_log) ? Number.NEGATIVE_INFINITY : 0.0; // R_D__0
417 }
418 if (!jstat.isFinite(lambda)) return (give_log) ? Number.NEGATIVE_INFINITY : 0.0; //R_D__0;
419 if (x < 0) return(give_log) ? Number.NEGATIVE_INFINITY : 0.0; //R_D__0
420 if (x <= lambda * jstat.DBL_MIN) {
421 return (give_log) ? -lambda : Math.exp(-lambda); // R_D_exp(-lambda)
422 }
423 if (lambda < x * jstat.DBL_MIN) {
424 var param = -lambda + x*Math.log(lambda) -jstat.lgamma(x+1);
425 return (give_log) ? param : Math.exp(param); // R_D_exp(-lambda + x*log(lambda) -lgammafn(x+1))
426 }
427 var param1 = jstat.TWO_PI * x; // f
428 var param2 = -jstat.stirlerr(x)-jstat.bd0(x,lambda); // x
429 return (give_log) ? -0.5*Math.log(param1)+param2 : Math.exp(param2)/Math.sqrt(param1); // R_D_fexp(M_2PI*x, -stirlerr(x)-bd0(x,lambda))
430 //return(R_D_fexp( , -stirlerr(x)-bd0(x,lambda) ));
431 }
432
433 /** Evaluates the "deviance part"
434 * bd0(x,M) := M * D0(x/M) = M*[ x/M * log(x/M) + 1 - (x/M) ] =
435 * = x * log(x/M) + M - x
436 * where M = E[X] = n*p (or = lambda), for x, M > 0
437 *
438 * in a manner that should be stable (with small relative error)
439 * for all x and M=np. In particular for x/np close to 1, direct
440 * evaluation fails, and evaluation is based on the Taylor series
441 * of log((1+v)/(1-v)) with v = (x-np)/(x+np).
442 */
443 jstat.bd0 = function(x, np) {
444 var ej, s, s1, v, j;
445 if(!jstat.isFinite(x) || !jstat.isFinite(np) || np == 0.0) throw "illegal parameter in jstat.bd0";
446
447 if(Math.abs(x-np) > 0.1*(x+np)) {
448 v = (x-np)/(x+np);
449 s = (x-np)*v;/* s using v -- change by MM */
450 ej = 2*x*v;
451 v = v*v;
452 for (j=1; ; j++) { /* Taylor series */
453 ej *= v;
454 s1 = s+ej/((j<<1)+1);
455 if (s1==s) /* last term was effectively 0 */
456 return(s1);
457 s = s1;
458 }
459 }
460 /* else: | x - np | is not too small */
461 return(x*Math.log(x/np)+np-x);
462 }
463
464 /** Computes the log of the error term in Stirling's formula.
465 * For n > 15, uses the series 1/12n - 1/360n^3 + ...
466 * For n <=15, integers or half-integers, uses stored values.
467 * For other n < 15, uses lgamma directly (don't use this to
468 * write lgamma!)
469 */
470 jstat.stirlerr= function(n) {
471 var S0 = 0.083333333333333333333;
472 var S1 = 0.00277777777777777777778;
473 var S2 = 0.00079365079365079365079365;
474 var S3 = 0.000595238095238095238095238;
475 var S4 = 0.0008417508417508417508417508;
476
477 var sferr_halves = [
478 0.0, /* n=0 - wrong, place holder only */
479 0.1534264097200273452913848, /* 0.5 */
480 0.0810614667953272582196702, /* 1.0 */
481 0.0548141210519176538961390, /* 1.5 */
482 0.0413406959554092940938221, /* 2.0 */
483 0.03316287351993628748511048, /* 2.5 */
484 0.02767792568499833914878929, /* 3.0 */
485 0.02374616365629749597132920, /* 3.5 */
486 0.02079067210376509311152277, /* 4.0 */
487 0.01848845053267318523077934, /* 4.5 */
488 0.01664469118982119216319487, /* 5.0 */
489 0.01513497322191737887351255, /* 5.5 */
490 0.01387612882307074799874573, /* 6.0 */
491 0.01281046524292022692424986, /* 6.5 */
492 0.01189670994589177009505572, /* 7.0 */
493 0.01110455975820691732662991, /* 7.5 */
494 0.010411265261972096497478567, /* 8.0 */
495 0.009799416126158803298389475, /* 8.5 */
496 0.009255462182712732917728637, /* 9.0 */
497 0.008768700134139385462952823, /* 9.5 */
498 0.008330563433362871256469318, /* 10.0 */
499 0.007934114564314020547248100, /* 10.5 */
500 0.007573675487951840794972024, /* 11.0 */
501 0.007244554301320383179543912, /* 11.5 */
502 0.006942840107209529865664152, /* 12.0 */
503 0.006665247032707682442354394, /* 12.5 */
504 0.006408994188004207068439631, /* 13.0 */
505 0.006171712263039457647532867, /* 13.5 */
506 0.005951370112758847735624416, /* 14.0 */
507 0.005746216513010115682023589, /* 14.5 */
508 0.005554733551962801371038690 /* 15.0 */
509 ];
510
511 var nn;
512
513 if (n <= 15.0) {
514 nn = n + n;
515 if (nn == parseInt(nn)) return(sferr_halves[parseInt(nn)]);
516 return(jstat.lgamma(n + 1.0) - (n + 0.5)*Math.log(n) + n - jstat.LN_SQRT_2PI);
517 }
518
519 nn = n*n;
520 if (n>500) return((S0-S1/nn)/n);
521 if (n> 80) return((S0-(S1-S2/nn)/nn)/n);
522 if (n> 35) return((S0-(S1-(S2-S3/nn)/nn)/nn)/n);
523 /* 15 < n <= 35 : */
524 return((S0-(S1-(S2-(S3-S4/nn)/nn)/nn)/nn)/n);
525 }
526
527
528
529 /** The function lgamma computes log|gamma(x)|. The function
530 * lgammafn_sign in addition assigns the sign of the gamma function
531 * to the address in the second argument if this is not null.
532 */
533 jstat.lgamma = function(x) {
534 function lgammafn_sign(x, sgn) {
535 var ans, y, sinpiy;
536 var xmax = 2.5327372760800758e+305;
537 var dxrel = 1.490116119384765696e-8;
538
539 // if (xmax == 0) {/* initialize machine dependent constants _ONCE_ */
540 // xmax = jstat.DBL_MAX/Math.log(jstat.DBL_MAX);/* = 2.533 e305 for IEEE double */
541 // dxrel = Math.sqrt(jstat.DBL_EPSILON);/* sqrt(Eps) ~ 1.49 e-8 for IEEE double */
542 // }
543
544 /* For IEEE double precision DBL_EPSILON = 2^-52 = 2.220446049250313e-16 :
545 xmax = DBL_MAX / log(DBL_MAX) = 2^1024 / (1024 * log(2)) = 2^1014 / log(2)
546 dxrel = sqrt(DBL_EPSILON) = 2^-26 = 5^26 * 1e-26 (is *exact* below !)
547 */
548
549 if (sgn != null) sgn = 1;
550
551 if(isNaN(x)) return x;
552
553 if (x < 0 && (Math.floor(-x) % 2.0) == 0)
554 if (sgn != null) sgn = -1;
555
556 if (x <= 0 && x == jstat.trunc(x)) { /* Negative integer argument */
557 console.warn("Negative integer argument in lgammafn_sign");
558 return Number.POSITIVE_INFINITY;/* +Inf, since lgamma(x) = log|gamma(x)| */
559 }
560
561 y = Math.abs(x);
562
563 if(y <= 10) return Math.log(Math.abs(jstat.gamma(x))); // TODO: implement jstat.gamma
564
565 if(y > xmax) {
566 console.warn("Illegal arguement passed to lgammafn_sign");
567 return Number.POSITIVE_INFINITY;
568 }
569
570 if(x > 0) {
571 if(x > 1e17) {
572 return (x*(Math.log(x)-1.0));
573 } else if(x > 4934720.0) {
574 return (jstat.LN_SQRT_2PI + (x-0.5) * Math.log(x) - x);
575 } else {
576 return jstat.LN_SQRT_2PI + (x-0.5) * Math.log(x) - x + jstat.lgammacor(x); // TODO: implement lgammacor
577 }
578 }
579
580 sinpiy = Math.abs(Math.sin(Math.PI * y));
581
582 if(sinpiy == 0) {
583 throw "Should never happen!!";
584 }
585
586 ans = jstat.LN_SQRT_PId2 + (x - 0.5) * Math.log(y) - x - Math.log(sinpiy) - jstat.lgammacor(y);
587
588 if(Math.abs((x-jstat.trunc(x-0.5))* ans / x) < dxrel) {
589 throw "The answer is less than half the precision argument too close to a negative integer";
590 }
591 return ans;
592 }
593
594 return lgammafn_sign(x, null);
595 }
596
597 jstat.gamma = function(x) {
598 var xbig = 171.624;
599 var p = [
600 -1.71618513886549492533811,
601 24.7656508055759199108314,-379.804256470945635097577,
602 629.331155312818442661052,866.966202790413211295064,
603 -31451.2729688483675254357,-36144.4134186911729807069,
604 66456.1438202405440627855
605 ];
606 var q = [
607 -30.8402300119738975254353,
608 315.350626979604161529144,-1015.15636749021914166146,
609 -3107.77167157231109440444,22538.1184209801510330112,
610 4755.84627752788110767815,-134659.959864969306392456,
611 -115132.259675553483497211
612 ];
613 var c = [
614 -.001910444077728,8.4171387781295e-4,
615 -5.952379913043012e-4,7.93650793500350248e-4,
616 -.002777777777777681622553,.08333333333333333331554247,
617 .0057083835261
618 ];
619
620 var i,n,parity,fact,xden,xnum,y,z,yi,res,sum,ysq;
621
622 parity = (0);
623 fact = 1.0;
624 n = 0;
625 y=x;
626 if(y <= 0.0) {
627 /* -------------------------------------------------------------
628 Argument is negative
629 ------------------------------------------------------------- */
630 y = -x;
631 yi = jstat.trunc(y);
632 res = y - yi;
633 if (res != 0.0) {
634 if (yi != jstat.trunc(yi * 0.5) * 2.0)
635 parity = (1);
636 fact = -Math.PI / Math.sin(Math.PI * res);
637 y += 1.0;
638 } else {
639 return(Number.POSITIVE_INFINITY);
640 }
641 }
642 /* -----------------------------------------------------------------
643 Argument is positive
644 -----------------------------------------------------------------*/
645 if (y < jstat.DBL_EPSILON) {
646 /* --------------------------------------------------------------
647 Argument < EPS
648 -------------------------------------------------------------- */
649 if (y >= jstat.DBL_MIN) {
650 res = 1.0 / y;
651 } else {
652 return(Number.POSITIVE_INFINITY);
653 }
654 } else if (y < 12.0) {
655 yi = y;
656 if (y < 1.0) {
657 /* ---------------------------------------------------------
658 EPS < argument < 1
659 --------------------------------------------------------- */
660 z = y;
661 y += 1.0;
662 } else {
663 /* -----------------------------------------------------------
664 1 <= argument < 12, reduce argument if necessary
665 ----------------------------------------------------------- */
666 n = parseInt(y) - 1;
667 y -= parseFloat(n);
668 z = y - 1.0;
669 }
670 /* ---------------------------------------------------------
671 Evaluate approximation for 1. < argument < 2.
672 ---------------------------------------------------------*/
673 xnum = 0.0;
674 xden = 1.0;
675 for (i = 0; i < 8; ++i) {
676 xnum = (xnum + p[i]) * z;
677 xden = xden * z + q[i];
678 }
679 res = xnum / xden + 1.0;
680 if (yi < y) {
681 /* --------------------------------------------------------
682 Adjust result for case 0. < argument < 1.
683 -------------------------------------------------------- */
684 res /= yi;
685 } else if (yi > y) {
686 /* ----------------------------------------------------------
687 Adjust result for case 2. < argument < 12.
688 ---------------------------------------------------------- */
689 for (i = 0; i < n; ++i) {
690 res *= y;
691 y += 1.0;
692 }
693 }
694 } else {
695 /* -------------------------------------------------------------
696 Evaluate for argument >= 12.,
697 ------------------------------------------------------------- */
698 if (y <= xbig) {
699 ysq = y * y;
700 sum = c[6];
701 for (i = 0; i < 6; ++i) {
702 sum = sum / ysq + c[i];
703 }
704 sum = sum / y - y + jstat.LN_SQRT_2PI;
705 sum += (y - 0.5) * Math.log(y);
706 res = Math.exp(sum);
707 } else {
708 return(Number.POSITIVE_INFINITY);
709 }
710 }
711 /* ----------------------------------------------------------------------
712 Final adjustments and return
713 ----------------------------------------------------------------------*/
714 if (parity)
715 res = -res;
716 if (fact != 1.0)
717 res = fact / res;
718 return res;
719 }
720
721 /** Compute the log gamma correction factor for x >= 10 so that
722 *
723 * log(gamma(x)) = .5*log(2*pi) + (x-.5)*log(x) -x + lgammacor(x)
724 *
725 * [ lgammacor(x) is called Del(x) in other contexts (e.g. dcdflib)]
726 */
727 jstat.lgammacor = function(x) {
728 var algmcs = [
729 +.1666389480451863247205729650822e+0,
730 -.1384948176067563840732986059135e-4,
731 +.9810825646924729426157171547487e-8,
732 -.1809129475572494194263306266719e-10,
733 +.6221098041892605227126015543416e-13,
734 -.3399615005417721944303330599666e-15,
735 +.2683181998482698748957538846666e-17,
736 -.2868042435334643284144622399999e-19,
737 +.3962837061046434803679306666666e-21,
738 -.6831888753985766870111999999999e-23,
739 +.1429227355942498147573333333333e-24,
740 -.3547598158101070547199999999999e-26,
741 +.1025680058010470912000000000000e-27,
742 -.3401102254316748799999999999999e-29,
743 +.1276642195630062933333333333333e-30
744 ];
745
746 var tmp;
747 var nalgm = 5;
748 var xbig = 94906265.62425156;
749 var xmax = 3.745194030963158e306;
750
751 if(x < 10) {
752 return Number.NaN;
753 } else if (x >= xmax) {
754 throw "Underflow error in lgammacor";
755 } else if (x < xbig) {
756 tmp = 10 / x;
757 return jstat.chebyshev(tmp*tmp*2-1,algmcs,nalgm) / x;
758 }
759 return 1 / (x*12);
760 }
761
762 /*
763 * Incomplete Beta function
764 */
765 jstat.incompleteBeta = function(a, b, x) {
766 /*
767 * Used by incompleteBeta: Evaluates continued fraction for incomplete
768 * beta function by modified Lentz's method.
769 */
770 function betacf(a, b, x) {
771 var MAXIT = 100;
772 var EPS = 3.0e-12;
773 var FPMIN = 1.0e-30;
774
775 var m,m2,aa,c,d,del,h,qab,qam,qap;
776
777 qab=a+b;
778 qap=a+1.0;
779 qam=a-1.0;
780 c=1.0;
781 d=1.0-qab*x/qap;
782
783 if(Math.abs(d) < FPMIN) {
784 d=FPMIN;
785 }
786
787 d = 1.0/d;
788 h=d;
789 for(m = 1; m <= MAXIT; m++) {
790 m2=2*m;
791 aa=m*(b-m)*x/((qam+m2)*(a+m2));
792 d=1.0+aa*d;
793 if(Math.abs(d) < FPMIN) {
794 d = FPMIN;
795 }
796 c=1.0+aa/c;
797 if(Math.abs(c) < FPMIN) {
798 c = FPMIN;
799 }
800 d=1.0/d;
801 h *= d*c;
802 aa = -(a+m)*(qab+m)*x/((a+m2) * (qap+m2));
803 d=1.0+aa*d;
804 if(Math.abs(d) < FPMIN) {
805 d = FPMIN;
806 }
807 c=1.0+aa/c;
808 if(Math.abs(c) < FPMIN) {
809 c=FPMIN;
810 }
811 d=1.0/d;
812 del=d*c;
813 h *= del;
814 if(Math.abs(del-1.0) < EPS) {
815 // are we done?
816 break;
817 }
818 }
819 if(m > MAXIT) {
820 console.warn("a or b too big, or MAXIT too small in betacf: " + a + ", " + b + ", " + x + ", " + h);
821 return h;
822 }
823 if(isNaN(h)) {
824 console.warn(a + ", " + b + ", " + x);
825 }
826 return h;
827 }
828
829 var bt;
830
831 if(x < 0.0 || x > 1.0) {
832 throw "bad x in routine incompleteBeta";
833 }
834 if(x == 0.0 || x == 1.0) {
835 bt = 0.0;
836 } else {
837 bt = Math.exp(jstat.lgamma(a+b) - jstat.lgamma(a) - jstat.lgamma(b) + a * Math.log(x)+ b * Math.log(1.0-x));
838 }
839 if(x < (a + 1.0)/(a+b+2.0)) {
840 return bt * betacf(a,b,x)/a;
841 } else {
842 return 1.0-bt*betacf(b,a,1.0-x)/b;
843 }
844 }
845
846 /** Evaluates the n-term Chebyshev series
847 * "a" at "x".
848 */
849 jstat.chebyshev = function(x, a, n) {
850 var b0, b1, b2, twox;
851 var i;
852
853 if (n < 1 || n > 1000) return Number.NaN;
854
855 if (x < -1.1 || x > 1.1) return Number.NaN;
856
857 twox = x * 2;
858 b2 = b1 = 0;
859 b0 = 0;
860 for (i = 1; i <= n; i++) {
861 b2 = b1;
862 b1 = b0;
863 b0 = twox * b1 - b2 + a[n - i];
864 }
865 return (b0 - b2) * 0.5;
866 }
867
868 jstat.fmin2 = function(x, y) {
869 return (x < y) ? x : y;
870 }
871
872 jstat.log1p = function(x) {
873 // http://kevin.vanzonneveld.net
874 // + original by: Brett Zamir (http://brett-zamir.me)
875 // % note 1: Precision 'n' can be adjusted as desired
876 // * example 1: log1p(1e-15);
877 // * returns 1: 9.999999999999995e-16
878
879 var ret = 0,
880 n = 50; // degree of precision
881 if (x <= -1) {
882 return Number.NEGATIVE_INFINITY; // JavaScript style would be to return Number.NEGATIVE_INFINITY
883 }
884 if (x < 0 || x > 1) {
885 return Math.log(1 + x);
886 }
887 for (var i = 1; i < n; i++) {
888 if ((i % 2) === 0) {
889 ret -= Math.pow(x, i) / i;
890 } else {
891 ret += Math.pow(x, i) / i;
892 }
893 }
894 return ret;
895 }
896
897 jstat.expm1 = function(x) {
898 var y, a = Math.abs(x);
899 if(a < jstat.DBL_EPSILON) return x;
900 if(a > 0.697) return Math.exp(x) - 1; /* negligable cancellation */
901 if(a > 1e-8) {
902 y = Math.exp(x) - 1;
903 } else {
904 y = (x / 2 + 1) * x;
905 }
906
907 /* Newton step for solving log(1 + y) = x for y : */
908 /* WARNING: does not work for y ~ -1: bug in 1.5.0 */
909 y -= (1 + y) * (jstat.log1p(y) - x);
910 return y;
911 }
912
913 jstat.logBeta = function(a, b) {
914 var corr, p, q;
915 p = q = a;
916 if(b < p) p = b;/* := min(a,b) */
917 if(b > q) q = b;/* := max(a,b) */
918
919 /* both arguments must be >= 0 */
920 if (p < 0) {
921 console.warn('Both arguements must be >= 0');
922 return Number.NaN;
923 }
924 else if (p == 0) {
925 return Number.POSITIVE_INFINITY;
926 }
927 else if (!jstat.isFinite(q)) { /* q == +Inf */
928 return Number.NEGATIVE_INFINITY;
929 }
930
931 if (p >= 10) {
932 /* p and q are big. */
933 corr = jstat.lgammacor(p) + jstat.lgammacor(q) - jstat.lgammacor(p + q);
934 return Math.log(q) * -0.5 + jstat.LN_SQRT_2PI + corr
935 + (p - 0.5) * Math.log(p / (p + q)) + q * jstat.log1p(-p / (p + q));
936 }
937 else if (q >= 10) {
938 /* p is small, but q is big. */
939 corr = jstat.lgammacor(q) - jstat.lgammacor(p + q);
940 return jstat.lgamma(p) + corr + p - p * Math.log(p + q)
941 + (q - 0.5) * jstat.log1p(-p / (p + q));
942 }
943 else
944 /* p and q are small: p <= q < 10. */
945 return Math.log(jstat.gamma(p) * (jstat.gamma(q) / jstat.gamma(p + q)));
946 }
947
948 jstat.dbinom_raw = function(x, n, p, q, give_log) {
949 if(give_log == null) give_log = false;
950 var lf, lc;
951
952 if(p == 0) {
953 if(x == 0) {
954 // R_D__1
955 return (give_log) ? 0.0 : 1.0;
956 } else {
957 // R_D__0
958 return (give_log) ? Number.NEGATIVE_INFINITY : 0.0;
959 }
960 }
961 if(q == 0) {
962 if(x == n) {
963 // R_D__1
964 return (give_log) ? 0.0 : 1.0;
965 } else {
966 // R_D__0
967 return (give_log) ? Number.NEGATIVE_INFINITY : 0.0;
968 }
969 }
970
971 if (x == 0) {
972 if(n == 0) return (give_log) ? 0.0 : 1.0; //R_D__1;
973 lc = (p < 0.1) ? -jstat.bd0(n,n*q) - n*p : n*Math.log(q);
974 return ( give_log ) ? lc : Math.exp(lc); //R_D_exp(lc)
975 }
976
977 if (x == n) {
978 lc = (q < 0.1) ? -jstat.bd0(n,n*p) - n*q : n*Math.log(p);
979 return ( give_log ) ? lc : Math.exp(lc); //R_D_exp(lc)
980 }
981
982 if (x < 0 || x > n) return (give_log) ? Number.NEGATIVE_INFINITY : 0.0; // R_D__0;
983
984 /* n*p or n*q can underflow to zero if n and p or q are small. This
985 used to occur in dbeta, and gives NaN as from R 2.3.0. */
986 lc = jstat.stirlerr(n) - jstat.stirlerr(x) - jstat.stirlerr(n-x) - jstat.bd0(x,n*p) - jstat.bd0(n-x,n*q);
987
988 /* f = (M_2PI*x*(n-x))/n; could overflow or underflow */
989 /* Upto R 2.7.1:
990 * lf = log(M_2PI) + log(x) + log(n-x) - log(n);
991 * -- following is much better for x << n : */
992 lf = Math.log(jstat.TWO_PI) + Math.log(x) + jstat.log1p(- x/n);
993
994 return (give_log) ? lc - 0.5*lf : Math.exp(lc - 0.5*lf); // R_D_exp(lc - 0.5*lf);
995 }
996
997 jstat.max = function(values) {
998 var max = Number.NEGATIVE_INFINITY;
999 for(var i = 0; i < values.length; i++) {
1000 if(values[i] > max) {
1001 max = values[i];
1002 }
1003 }
1004 return max;
1005 }
1006
1007 /******************************************************************************/
1008 /* Probability Distributions */
1009 /******************************************************************************/
1010
1011 /**
1012 * Range class
1013 */
1014 var Range = Class.extend({
1015 init: function(min, max, numPoints) {
1016 this._minimum = parseFloat(min);
1017 this._maximum = parseFloat(max);
1018 this._numPoints = parseFloat(numPoints);
1019 },
1020 getMinimum: function() {
1021 return this._minimum;
1022 },
1023 getMaximum: function() {
1024 return this._maximum;
1025 },
1026 getNumPoints: function() {
1027 return this._numPoints;
1028 },
1029 getPoints: function() {
1030 var results = [];
1031 var x = this._minimum;
1032 var step = (this._maximum-this._minimum)/(this._numPoints-1);
1033 for(var i = 0; i < this._numPoints; i++) {
1034 results[i] = parseFloat(x.toFixed(6));
1035 x += step;
1036 }
1037 return results;
1038 }
1039 });
1040
1041 Range.validate = function(range) {
1042 if( ! range instanceof Range) {
1043 return false;
1044 }
1045 if(isNaN(range.getMinimum()) || isNaN(range.getMaximum()) || isNaN(range.getNumPoints()) || range.getMaximum() < range.getMinimum() || range.getNumPoints() <= 0) {
1046 return false;
1047 }
1048 return true;
1049 }
1050
1051 var ContinuousDistribution = Class.extend({
1052 init: function(name) {
1053 this._name = name;
1054 },
1055 toString: function() {
1056 return this._string;
1057 },
1058 getName: function() {
1059 return this._name;
1060 },
1061 getClassName: function() {
1062 return this._name + 'Distribution';
1063 },
1064 density: function(valueOrRange) {
1065 if(!isNaN(valueOrRange)) {
1066 // single value
1067 return parseFloat(this._pdf(valueOrRange).toFixed(15));
1068 } else if (Range.validate(valueOrRange)) {
1069 // multiple values
1070 var points = valueOrRange.getPoints();
1071
1072 var result = [];
1073 // For each point in the range
1074 for(var i = 0; i < points.length; i++) {
1075 result[i] = parseFloat(this._pdf(points[i]));
1076 }
1077 return result;
1078 } else {
1079 // neither value or range
1080 throw "Invalid parameter supplied to " + this.getClassName() + ".density()";
1081 }
1082 },
1083 cumulativeDensity: function(valueOrRange) {
1084 if(!isNaN(valueOrRange)) {
1085 // single value
1086 return parseFloat(this._cdf(valueOrRange).toFixed(15));
1087 } else if (Range.validate(valueOrRange)) {
1088 // multiple values
1089 var points = valueOrRange.getPoints();
1090 var result = [];
1091 // For each point in the range
1092 for(var i = 0; i < points.length; i++) {
1093 result[i] = parseFloat(this._cdf(points[i]));
1094 }
1095 return result;
1096 } else {
1097 // neither value or range
1098 throw "Invalid parameter supplied to " + this.getClassName() + ".cumulativeDensity()";
1099 }
1100 },
1101 getRange: function(standardDeviations, numPoints) {
1102 if(standardDeviations == null) {
1103 standardDeviations = 5;
1104 }
1105 if(numPoints == null) {
1106 numPoints = 100;
1107 }
1108 var min = this.getMean() - standardDeviations * Math.sqrt(this.getVariance());
1109 var max = this.getMean() + standardDeviations * Math.sqrt(this.getVariance());
1110
1111 if(this.getClassName() == 'GammaDistribution' || this.getClassName() == 'LogNormalDistribution') {
1112 min = 0.0;
1113 max = this.getMean() + standardDeviations * Math.sqrt(this.getVariance());
1114 } else if(this.getClassName() == 'BetaDistribution') {
1115 min = 0.0;
1116 max = 1.0;
1117 }
1118
1119
1120 var range = new Range(min, max, numPoints);
1121 return range;
1122 },
1123 getVariance: function(){},
1124 getMean: function(){},
1125 getQuantile: function(p) {
1126 var self = this;
1127 /*
1128 * Recursive function to find the closest match
1129 */
1130 function findClosestMatch(range, p) {
1131 var ERR = 1.0e-5;
1132 var xs = range.getPoints();
1133 var closestIndex = 0;
1134 var closestDistance = 999;
1135
1136 for(var i=0; i<xs.length; i++) {
1137 var pp = self.cumulativeDensity(xs[i]);
1138 var distance = Math.abs(pp - p);
1139 if(distance < closestDistance) {
1140 // closer value found
1141 closestIndex = i;
1142 closestDistance = distance;
1143 }
1144 }
1145 if(closestDistance <= ERR) {
1146 // Acceptable - return value;
1147 return xs[closestIndex];
1148 } else {
1149 // Calculate the new range
1150 var newRange = new Range(xs[closestIndex-1], xs[closestIndex+1],20);
1151 return findClosestMatch(newRange, p);
1152 }
1153 }
1154 var range = this.getRange(5, 20);
1155 return findClosestMatch(range, p);
1156 }
1157 });
1158
1159 /**
1160 * A normal distribution object
1161 */
1162 var NormalDistribution = ContinuousDistribution.extend({
1163 init: function(mean, sigma) {
1164 this._super('Normal');
1165 this._mean = parseFloat(mean);
1166 this._sigma = parseFloat(sigma);
1167 this._string = "Normal ("+this._mean.toFixed(2)+", " + this._sigma.toFixed(2) + ")";
1168 },
1169 _pdf: function(x, give_log) {
1170 if(give_log == null) {
1171 give_log=false;
1172 } // default is false;
1173 var sigma = this._sigma;
1174 var mu = this._mean;
1175 if(!jstat.isFinite(sigma)) {
1176 return (give_log) ? Number.NEGATIVE_INFINITY : 0.0
1177 }
1178 if(!jstat.isFinite(x) && mu == x) {
1179 return Number.NaN;
1180 }
1181 if(sigma<=0) {
1182 if(sigma < 0) {
1183 throw "invalid sigma in _pdf";
1184 }
1185 return (x==mu)?Number.POSITIVE_INFINITY:(give_log)?Number.NEGATIVE_INFINITY:0.0;
1186 }
1187 x=(x-mu)/sigma;
1188 if(!jstat.isFinite(x)){
1189 return (give_log)?Number.NEGATIVE_INFINITY:0.0;
1190 }
1191 return (give_log ? -(jstat.LN_SQRT_2PI + 0.5 * x * x + Math.log(sigma)) :
1192 jstat.ONE_SQRT_2PI * Math.exp(-0.5 * x * x) / sigma);
1193 },
1194 _cdf: function(x, lower_tail, log_p) {
1195
1196 if(lower_tail == null) lower_tail = true;
1197 if(log_p == null) log_p = false;
1198
1199 function pnorm_both(x, cum, ccum, i_tail, log_p) {
1200 /* i_tail in {0,1,2} means: "lower", "upper", or "both" :
1201 if(lower) return *cum := P[X <= x]
1202 if(upper) return *ccum := P[X > x] = 1 - P[X <= x]
1203 */
1204
1205 var a = [
1206 2.2352520354606839287,
1207 161.02823106855587881,
1208 1067.6894854603709582,
1209 18154.981253343561249,
1210 0.065682337918207449113
1211 ];
1212 var b = [
1213 47.20258190468824187,
1214 976.09855173777669322,
1215 10260.932208618978205,
1216 45507.789335026729956
1217 ];
1218 var c = [
1219 0.39894151208813466764,
1220 8.8831497943883759412,
1221 93.506656132177855979,
1222 597.27027639480026226,
1223 2494.5375852903726711,
1224 6848.1904505362823326,
1225 11602.651437647350124,
1226 9842.7148383839780218,
1227 1.0765576773720192317e-8
1228 ];
1229 var d = [
1230 22.266688044328115691,
1231 235.38790178262499861,
1232 1519.377599407554805,
1233 6485.558298266760755,
1234 18615.571640885098091,
1235 34900.952721145977266,
1236 38912.003286093271411,
1237 19685.429676859990727
1238 ];
1239 var p = [
1240 0.21589853405795699,
1241 0.1274011611602473639,
1242 0.022235277870649807,
1243 0.001421619193227893466,
1244 2.9112874951168792e-5,
1245 0.02307344176494017303
1246 ];
1247 var q = [
1248 1.28426009614491121,
1249 0.468238212480865118,
1250 0.0659881378689285515,
1251 0.00378239633202758244,
1252 7.29751555083966205e-5
1253 ];
1254
1255 var xden, xnum, temp, del, eps, xsq, y, i, lower, upper;
1256
1257 /* Consider changing these : */
1258 eps = jstat.DBL_EPSILON * 0.5;
1259
1260 /* i_tail in {0,1,2} =^= {lower, upper, both} */
1261 lower = i_tail != 1;
1262 upper = i_tail != 0;
1263
1264 y = Math.abs(x);
1265
1266 if (y <= 0.67448975) { /* qnorm(3/4) = .6744.... -- earlier had 0.66291 */
1267 if (y > eps) {
1268 xsq = x * x;
1269 xnum = a[4] * xsq;
1270 xden = xsq;
1271 for (i = 0; i < 3; ++i) {
1272 xnum = (xnum + a[i]) * xsq;
1273 xden = (xden + b[i]) * xsq;
1274 }
1275 } else {
1276 xnum = xden = 0.0;
1277 }
1278 temp = x * (xnum + a[3]) / (xden + b[3]);
1279 if(lower) cum = 0.5 + temp;
1280 if(upper) ccum = 0.5 - temp;
1281 if(log_p) {
1282 if(lower) cum = Math.log(cum);
1283 if(upper) ccum = Math.log(ccum);
1284 }
1285
1286 } else if (y <= jstat.SQRT_32) {
1287 /* Evaluate pnorm for 0.674.. = qnorm(3/4) < |x| <= sqrt(32) ~= 5.657 */
1288
1289 xnum = c[8] * y;
1290 xden = y;
1291 for (i = 0; i < 7; ++i) {
1292 xnum = (xnum + c[i]) * y;
1293 xden = (xden + d[i]) * y;
1294 }
1295 temp = (xnum + c[7]) / (xden + d[7]);
1296
1297 /* do_del */
1298 xsq = jstat.trunc(x * 16) / 16;
1299 del = (x - xsq) * (x + xsq);
1300 if(log_p) {
1301 cum = (-xsq * xsq * 0.5) + (-del * 0.5) + Math.log(temp);
1302 if((lower && x > 0.) || (upper && x <= 0.))
1303 ccum = jstat.log1p(-Math.exp(-xsq * xsq * 0.5) *
1304 Math.exp(-del * 0.5) * temp);
1305 }
1306 else {
1307 cum = Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp;
1308 ccum = 1.0 - cum;
1309 }
1310 /* end do_del */
1311
1312 /* swap_tail */
1313 if (x > 0.0) {/* swap ccum <--> cum */
1314 temp = cum;
1315 if(lower) {
1316 cum = ccum;
1317
1318 }
1319 ccum = temp;
1320 }
1321 /* end swap_tail */
1322
1323 }
1324 /* else |x| > sqrt(32) = 5.657 :
1325 * the next two case differentiations were really for lower=T, log=F
1326 * Particularly *not* for log_p !
1327
1328 * Cody had (-37.5193 < x && x < 8.2924) ; R originally had y < 50
1329 *
1330 * Note that we do want symmetry(0), lower/upper -> hence use y
1331 */
1332
1333 else if((log_p && y < 1e170)|| (lower && -37.5193 < x && x < 8.2924)
1334 || (upper && -8.2924 < x && x < 37.5193)) {
1335 /* Evaluate pnorm for x in (-37.5, -5.657) union (5.657, 37.5) */
1336 xsq = 1.0 / (x * x); /* (1./x)*(1./x) might be better */
1337 xnum = p[5] * xsq;
1338 xden = xsq;
1339 for (i = 0; i < 4; ++i) {
1340 xnum = (xnum + p[i]) * xsq;
1341 xden = (xden + q[i]) * xsq;
1342 }
1343 temp = xsq * (xnum + p[4]) / (xden + q[4]);
1344 temp = (jstat.ONE_SQRT_2PI - temp) / y;
1345
1346 /* do_del */
1347 xsq = jstat.trunc(x * 16) / 16;
1348 del = (x - xsq) * (x + xsq);
1349 if(log_p) {
1350 cum = (-xsq * xsq * 0.5) + (-del * 0.5) + Math.log(temp);
1351 if((lower && x > 0.) || (upper && x <= 0.))
1352 ccum = jstat.log1p(-Math.exp(-xsq * xsq * 0.5) *
1353 Math.exp(-del * 0.5) * temp);
1354 }
1355 else {
1356 cum = Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp;
1357 ccum = 1.0 - cum;
1358 }
1359 /* end do_del */
1360
1361 /* swap_tail */
1362 if (x > 0.0) {/* swap ccum <--> cum */
1363 temp = cum;
1364 if(lower) {
1365 cum = ccum;
1366
1367 }
1368 ccum = temp;
1369 }
1370 /* end swap_tail */
1371
1372 } else { /* large x such that probs are 0 or 1 */
1373 if(x > 0) {
1374 cum = (log_p) ? 0.0 : 1.0; // R_D__1
1375 ccum = (log_p) ? Number.NEGATIVE_INFINITY : 0.0; //R_D__0;
1376 } else {
1377 cum = (log_p) ? Number.NEGATIVE_INFINITY : 0.0; //R_D__0;
1378 ccum = (log_p) ? 0.0 : 1.0; // R_D__1
1379 }
1380 }
1381
1382 return [cum, ccum];
1383 }
1384
1385 var p, cp;
1386 var mu = this._mean;
1387 var sigma = this._sigma;
1388 var R_DT_0, R_DT_1;
1389
1390 if(lower_tail) {
1391 if(log_p) {
1392 R_DT_0 = Number.NEGATIVE_INFINITY;
1393 R_DT_1 = 0.0;
1394 } else {
1395 R_DT_0 = 0.0;
1396 R_DT_1 = 1.0;
1397 }
1398 } else {
1399 if(log_p) {
1400 R_DT_0 = 0.0;
1401 R_DT_1 = Number.NEGATIVE_INFINITY;
1402 } else {
1403 R_DT_0 = 1.0;
1404 R_DT_1 = 0.0;
1405 }
1406 }
1407
1408 if(!jstat.isFinite(x) && mu == x) return Number.NaN;
1409 if(sigma <= 0) {
1410 if(sigma < 0) {
1411 console.warn("Sigma is less than 0");
1412 return Number.NaN;
1413 }
1414 return (x < mu) ? R_DT_0 : R_DT_1;
1415 }
1416
1417 p = (x - mu) / sigma;
1418
1419 if(!jstat.isFinite(p)) {
1420 return (x < mu) ? R_DT_0 : R_DT_1;
1421 }
1422
1423 x = p;
1424
1425 // pnorm_both(x, &p, &cp, (lower_tail ? 0 : 1), log_p);
1426 // result[0] == &p
1427 // result[1] == &cp
1428
1429 var result = pnorm_both(x, p, cp, (lower_tail ? false : true), log_p);
1430
1431 return (lower_tail ? result[0] : result[1]);
1432
1433 },
1434 getMean: function() {
1435 return this._mean;
1436 },
1437 getSigma: function() {
1438 return this._sigma;
1439 },
1440 getVariance: function() {
1441 return this._sigma*this._sigma;
1442 }
1443 });
1444
1445 /**
1446 * A Log-normal distribution object
1447 */
1448 var LogNormalDistribution = ContinuousDistribution.extend({
1449 init: function(location, scale) {
1450 this._super('LogNormal')
1451 this._location = parseFloat(location);
1452 this._scale = parseFloat(scale);
1453 this._string = "LogNormal ("+this._location.toFixed(2)+", " + this._scale.toFixed(2) + ")";
1454 },
1455 _pdf: function(x, give_log) {
1456 var y;
1457 var sdlog = this._scale;
1458 var meanlog = this._location;
1459 if(give_log == null) {
1460 give_log = false;
1461 }
1462
1463 if(sdlog <= 0) throw "Illegal parameter in _pdf";
1464
1465 if(x <= 0) {
1466 return (give_log) ? Number.NEGATIVE_INFINITY : 0.0;
1467 }
1468
1469 y = (Math.log(x) - meanlog) / sdlog;
1470
1471 return (give_log ? -(jstat.LN_SQRT_2PI + 0.5 * y * y + Math.log(x * sdlog)) :
1472 jstat.ONE_SQRT_2PI * Math.exp(-0.5 * y * y) / (x * sdlog));
1473
1474 },
1475 _cdf: function(x, lower_tail, log_p) {
1476 var sdlog = this._scale;
1477 var meanlog = this._location;
1478 if(lower_tail == null) {
1479 lower_tail = true;
1480 }
1481 if(log_p == null) {
1482 log_p = false;
1483 }
1484
1485
1486 if(sdlog <= 0) {
1487 throw "illegal std in _cdf";
1488 }
1489
1490 if(x > 0) {
1491 var nd = new NormalDistribution(meanlog, sdlog);
1492 return nd._cdf(Math.log(x), lower_tail, log_p);
1493 }
1494 if(lower_tail) {
1495 return (log_p) ? Number.NEGATIVE_INFINITY : 0.0; // R_D__0
1496 } else {
1497 return (log_p) ? 0.0 : 1.0; // R_D__1
1498 }
1499 },
1500 getLocation: function() {
1501 return this._location;
1502 },
1503 getScale: function() {
1504 return this._scale;
1505 },
1506 getMean: function() {
1507 return Math.exp((this._location + this._scale) / 2);
1508 },
1509 getVariance: function() {
1510 var ans = (Math.exp(this._scale)-1)*Math.exp(2*this._location+this._scale);
1511 return ans;
1512 }
1513 });
1514
1515
1516 /**
1517 * Gamma distribution object
1518 */
1519 var GammaDistribution = ContinuousDistribution.extend({
1520 init: function(shape, scale) {
1521 this._super('Gamma');
1522 this._shape = parseFloat(shape);
1523 this._scale = parseFloat(scale);
1524 this._string = "Gamma ("+this._shape.toFixed(2)+", " + this._scale.toFixed(2) + ")";
1525 },
1526 _pdf: function(x, give_log) {
1527 var pr;
1528 var shape = this._shape;
1529 var scale = this._scale;
1530 if(give_log == null) {
1531 give_log = false; // default value
1532 }
1533
1534 if(shape < 0 || scale <= 0) {
1535 throw "Illegal argument in _pdf";
1536 }
1537
1538 if(x < 0) {
1539 return (give_log) ? Number.NEGATIVE_INFINITY : 0.0; // R_D__0
1540 }
1541 if(shape == 0) { /* point mass at 0 */
1542 return (x == 0) ? Number.POSITIVE_INFINITY : (give_log) ? Number.NEGATIVE_INFINITY : 0.0; // R_D__0
1543 }
1544 if(x == 0) {
1545 if(shape < 1) return Number.POSITIVE_INFINITY;
1546 if(shape > 1) return (give_log) ? Number.NEGATIVE_INFINITY : 0.0; // R_D__0
1547 /* else */
1548 return (give_log) ? -Math.log(scale) : 1/scale;
1549 }
1550
1551 if(shape < 1) {
1552 pr = jstat.dopois_raw(shape, x/scale, give_log);
1553 return give_log ? pr + Math.log(shape/x) : pr*shape/x;
1554 }
1555 /* else shape >= 1 */
1556 pr = jstat.dopois_raw(shape-1, x/scale, give_log);
1557 return give_log ? pr - Math.log(scale) : pr/scale;
1558
1559 },
1560 /**
1561 * This function computes the distribution function for the
1562 * gamma distribution with shape parameter alph and scale parameter
1563 * scale. This is also known as the incomplete gamma function.
1564 * See Abramowitz and Stegun (6.5.1) for example.
1565 */
1566 _cdf: function(x, lower_tail, log_p) {
1567 /* define USE_PNORM */
1568 function USE_PNORM() {
1569 pn1 = Math.sqrt(alph) * 3.0 * (Math.pow(x/alph,1.0/3.0) + 1.0 / (9.0 * alph) - 1.0);
1570 var norm_dist = new NormalDistribution(0.0, 1.0);
1571 return norm_dist._cdf(pn1, lower_tail, log_p);
1572 }
1573
1574 /* Defaults */
1575 if(lower_tail == null) lower_tail = true;
1576 if(log_p == null) log_p = false;
1577 var alph = this._shape;
1578 var scale = this._scale;
1579 var xbig = 1.0e+8;
1580 var xlarge = 1.0e+37;
1581 var alphlimit = 1e5;
1582 var pn1,pn2,pn3,pn4,pn5,pn6,arg,a,b,c,an,osum,sum,n,pearson;
1583
1584 if(alph <= 0. || scale <= 0.) {
1585 console.warn('Invalid gamma params in _cdf');
1586 return Number.NaN;
1587 }
1588
1589 x/=scale;
1590 if(isNaN(x)) return x;
1591 if(x <= 0.0) {
1592 // R_DT_0
1593 if(lower_tail) {
1594 // R_D__0
1595 return (log_p) ? Number.NEGATIVE_INFINITY : 0.0;
1596 } else {
1597 // R_D__1
1598 return (log_p) ? 0.0 : 1.0;
1599 }
1600 }
1601
1602 if(alph > alphlimit) {
1603 return USE_PNORM();
1604 }
1605
1606 if(x > xbig * alph) {
1607 if(x > jstat.DBL_MAX * alph) {
1608 // R_DT_1
1609 if(lower_tail) {
1610 // R_D__1
1611 return (log_p) ? 0.0 : 1.0;
1612 } else {
1613 // R_D__0
1614 return (log_p) ? Number.NEGATIVE_INFINITY : 0.0;
1615 }
1616 } else {
1617 return USE_PNORM();
1618 }
1619 }
1620
1621 if(x <= 1.0 || x < alph) {
1622 pearson = 1; /* use pearson's series expansion */
1623 arg = alph * Math.log(x) - x - jstat.lgamma(alph + 1.0);
1624
1625 c = 1.0;
1626 sum = 1.0;
1627 a = alph;
1628 do {
1629 a += 1.0;
1630 c *= x / a;
1631 sum += c;
1632 } while(c > jstat.DBL_EPSILON * sum);
1633 } else { /* x >= max( 1, alph) */
1634 pearson = 0;/* use a continued fraction expansion */
1635 arg = alph * Math.log(x) - x - jstat.lgamma(alph);
1636
1637 a = 1. - alph;
1638 b = a + x + 1.;
1639 pn1 = 1.;
1640 pn2 = x;
1641 pn3 = x + 1.;
1642 pn4 = x * b;
1643 sum = pn3 / pn4;
1644
1645 for (n = 1; ; n++) {
1646 a += 1.;/* = n+1 -alph */
1647 b += 2.;/* = 2(n+1)-alph+x */
1648 an = a * n;
1649 pn5 = b * pn3 - an * pn1;
1650 pn6 = b * pn4 - an * pn2;
1651 if (Math.abs(pn6) > 0.) {
1652 osum = sum;
1653 sum = pn5 / pn6;
1654 if (Math.abs(osum - sum) <= jstat.DBL_EPSILON * jstat.fmin2(1.0, sum))
1655 break;
1656 }
1657 pn1 = pn3;
1658 pn2 = pn4;
1659 pn3 = pn5;
1660 pn4 = pn6;
1661
1662 if (Math.abs(pn5) >= xlarge) {
1663 pn1 /= xlarge;
1664 pn2 /= xlarge;
1665 pn3 /= xlarge;
1666 pn4 /= xlarge;
1667 }
1668 }
1669 }
1670 arg += Math.log(sum);
1671 lower_tail = (lower_tail == pearson);
1672
1673 if (log_p && lower_tail)
1674 return(arg);
1675 /* else */
1676 /* sum = exp(arg); and return if(lower_tail) sum else 1-sum : */
1677
1678 if(lower_tail) {
1679 return Math.exp(arg);
1680 } else {
1681 if(log_p) {
1682 // R_Log1_Exp(arg);
1683 return (arg > -Math.LN2) ? Math.log(-jstat.expm1(arg)) : jstat.log1p(-Math.exp(arg));
1684 } else {
1685 return -jstat.expm1(arg);
1686 }
1687 }
1688 },
1689 getShape: function() {
1690 return this._shape;
1691 },
1692 getScale: function() {
1693 return this._scale;
1694 },
1695 getMean: function() {
1696 return this._shape * this._scale;
1697 },
1698 getVariance: function() {
1699 return this._shape*Math.pow(this._scale,2);
1700 }
1701 });
1702
1703 /**
1704 * A Beta distribution object
1705 */
1706 var BetaDistribution = ContinuousDistribution.extend({
1707 init: function(alpha, beta) {
1708 this._super('Beta');
1709 this._alpha = parseFloat(alpha);
1710 this._beta = parseFloat(beta);
1711 this._string = "Beta ("+this._alpha.toFixed(2)+", " + this._beta.toFixed(2) + ")";
1712 },
1713 _pdf: function(x, give_log) {
1714 if(give_log == null) give_log = false; // default;
1715 var a = this._alpha;
1716 var b = this._beta;
1717 var lval;
1718 if(a <= 0 || b <= 0) {
1719 console.warn('Illegal arguments in _pdf');
1720 return Number.NaN;
1721 }
1722 if(x < 0 || x > 1) {
1723 // R_D__0
1724 return (give_log) ? Number.NEGATIVE_INFINITY : 0.0;
1725 }
1726 if(x == 0) {
1727 if(a > 1) {
1728 // R_D__0
1729 return (give_log) ? Number.NEGATIVE_INFINITY : 0.0;
1730 }
1731 if(a < 1) {
1732 return Number.POSITIVE_INFINITY;
1733 }
1734 /*a == 1 */ return (give_log) ? Math.log(b) : b; // R_D_val(b)
1735 }
1736 if(x == 1) {
1737 if(b > 1) {
1738 // R_D__0
1739 return (give_log) ? Number.NEGATIVE_INFINITY : 0.0;
1740 }
1741 if(b < 1) {
1742 return Number.POSITIVE_INFINITY;
1743 }
1744 /* b == 1 */ return (give_log) ? Math.log(a) : a; // R_D_val(a)
1745 }
1746 if(a<=2||b<=2) {
1747 lval = (a-1)*Math.log(x) + (b-1)*jstat.log1p(-x) - jstat.logBeta(a, b);
1748 } else {
1749 lval = Math.log(a+b-1) + jstat.dbinom_raw(a-1, a+b-2, x, 1-x, true);
1750 }
1751 //R_D_exp(lval)
1752 return (give_log) ? lval : Math.exp(lval);
1753 },
1754 _cdf: function(x, lower_tail, log_p) {
1755 if(lower_tail == null) lower_tail = true;
1756 if(log_p == null) log_p = false;
1757
1758 var pin = this._alpha;
1759 var qin = this._beta;
1760
1761 if(pin <= 0 || qin <= 0) {
1762 console.warn('Invalid argument in _cdf');
1763 return Number.NaN;
1764 }
1765
1766 if(x <= 0) {
1767 //R_DT_0;
1768 if(lower_tail) {
1769 // R_D__0
1770 return (log_p) ? Number.NEGATIVE_INFINITY : 0.0;
1771 } else {
1772 // R_D__1
1773 return (log_p) ? 0.1 : 1.0;
1774 }
1775 }
1776
1777 if(x >= 1){
1778 // R_DT_1
1779 if(lower_tail) {
1780 // R_D__1
1781 return (log_p) ? 0.1 : 1.0;
1782 } else {
1783 // R_D__0
1784 return (log_p) ? Number.NEGATIVE_INFINITY : 0.0;
1785 }
1786 }
1787
1788 /* else */
1789 return jstat.incompleteBeta(pin, qin, x);
1790 },
1791 getAlpha: function() {
1792 return this._alpha;
1793 },
1794 getBeta: function() {
1795 return this._beta;
1796 },
1797 getMean: function() {
1798 return this._alpha / (this._alpha+ this._beta);
1799 },
1800 getVariance: function() {
1801 var ans = (this._alpha * this._beta) / (Math.pow(this._alpha+this._beta,2)*(this._alpha+this._beta+1));
1802 return ans;
1803 }
1804 });
1805
1806 var StudentTDistribution = ContinuousDistribution.extend({
1807 init: function(degreesOfFreedom, mu) {
1808 this._super('StudentT');
1809 this._dof = parseFloat(degreesOfFreedom);
1810
1811 if(mu != null) {
1812 this._mu = parseFloat(mu);
1813 this._string = "StudentT ("+this._dof.toFixed(2)+", " + this._mu.toFixed(2)+ ")";
1814 } else {
1815 this._mu = 0.0;
1816 this._string = "StudentT ("+this._dof.toFixed(2)+")";
1817 }
1818
1819 },
1820 _pdf: function(x, give_log) {
1821 if(give_log == null) give_log = false;
1822 if(this._mu == null) {
1823 return this._dt(x, give_log);
1824 } else {
1825 var y = this._dnt(x, give_log);
1826 if(y > 1){
1827 console.warn('x:' + x + ', y: ' + y);
1828 }
1829 return y;
1830 }
1831 },
1832 _cdf: function(x, lower_tail, give_log) {
1833 if(lower_tail == null) lower_tail = true;
1834 if(give_log == null) give_log = false;
1835 if(this._mu == null) {
1836 return this._pt(x, lower_tail, give_log);
1837 } else {
1838 return this._pnt(x, lower_tail, give_log);
1839 }
1840 },
1841 _dt: function(x, give_log) {
1842 var t,u;
1843 var n = this._dof;
1844 if (n <= 0){
1845 console.warn('Invalid parameters in _dt');
1846 return Number.NaN;
1847 }
1848 if(!jstat.isFinite(x)) {
1849 return (give_log) ? Number.NEGATIVE_INFINITY : 0.0; // R_D__0;
1850 }
1851
1852 if(!jstat.isFinite(n)) {
1853 var norm = new NormalDistribution(0.0, 1.0);
1854 return norm.density(x, give_log);
1855 }
1856
1857
1858 t = -jstat.bd0(n/2.0,(n+1)/2.0) + jstat.stirlerr((n+1)/2.0) - jstat.stirlerr(n/2.0);
1859 if ( x*x > 0.2*n )
1860 u = Math.log( 1+ x*x/n ) * n/2;
1861 else
1862 u = -jstat.bd0(n/2.0,(n+x*x)/2.0) + x*x/2.0;
1863
1864 var p1 = jstat.TWO_PI *(1+x*x/n);
1865 var p2 = t-u;
1866
1867 return (give_log) ? -0.5*Math.log(p1) + p2 : Math.exp(p2)/Math.sqrt(p1); // R_D_fexp(M_2PI*(1+x*x/n), t-u);
1868 },
1869 _dnt: function(x, give_log) {
1870 if(give_log == null) give_log = false;
1871 var df = this._dof;
1872 var ncp = this._mu;
1873 var u;
1874
1875 if(df <= 0.0) {
1876 console.warn("Illegal arguments _dnf");
1877 return Number.NaN;
1878 }
1879 if(ncp == 0.0) {
1880 return this._dt(x, give_log);
1881 }
1882
1883 if(!jstat.isFinite(x)) {
1884 // R_D__0
1885 if(give_log) {
1886 return Number.NEGATIVE_INFINITY;
1887 } else {
1888 return 0.0;
1889 }
1890 }
1891
1892 /* If infinite df then the density is identical to a
1893 normal distribution with mean = ncp. However, the formula
1894 loses a lot of accuracy around df=1e9
1895 */
1896 if(!isFinite(df) || df > 1e8) {
1897 var dist = new NormalDistribution(ncp, 1.);
1898 return dist.density(x, give_log);
1899 }
1900
1901 /* Do calculations on log scale to stabilize */
1902
1903 /* Consider two cases: x ~= 0 or not */
1904 if (Math.abs(x) > Math.sqrt(df * jstat.DBL_EPSILON)) {
1905 var newT = new StudentTDistribution(df+2, ncp);
1906 u = Math.log(df) - Math.log(Math.abs(x)) +
1907 Math.log(Math.abs(newT._pnt(x*Math.sqrt((df+2)/df), true, false) -
1908 this._pnt(x, true, false)));
1909 /* FIXME: the above still suffers from cancellation (but not horribly) */
1910 }
1911 else { /* x ~= 0 : -> same value as for x = 0 */
1912 u = jstat.lgamma((df+1)/2) - jstat.lgamma(df/2)
1913 - .5*(Math.log(Math.PI) + Math.log(df) + ncp*ncp);
1914 }
1915
1916 return (give_log ? u : Math.exp(u));
1917 },
1918 _pt: function(x, lower_tail, log_p) {
1919 if(lower_tail == null) lower_tail = true;
1920 if(log_p == null) log_p = false;
1921 var val, nx;
1922 var n = this._dof;
1923 var DT_0, DT_1;
1924
1925 if(lower_tail) {
1926 if(log_p) {
1927 DT_0 = Number.NEGATIVE_INFINITY;
1928 DT_1 = 1.;
1929 } else {
1930 DT_0 = 0.;
1931 DT_1 = 1.;
1932 }
1933 } else {
1934 if(log_p) {
1935 // not lower_tail but log_p
1936 DT_0 = 0.;
1937 DT_1 = Number.NEGATIVE_INFINITY;
1938 } else {
1939 // not lower_tail and not log_p
1940 DT_0 = 1.;
1941 DT_1 = 0.;
1942 }
1943 }
1944
1945 if(n <= 0.0) {
1946 console.warn("Invalid T distribution _pt");
1947 return Number.NaN;
1948 }
1949 var norm = new NormalDistribution(0,1);
1950 if(!jstat.isFinite(x)) {
1951 return (x < 0) ? DT_0 : DT_1;
1952 }
1953 if(!jstat.isFinite(n)) {
1954 return norm._cdf(x, lower_tail, log_p);
1955 }
1956
1957 if (n > 4e5) { /*-- Fixme(?): test should depend on `n' AND `x' ! */
1958 /* Approx. from Abramowitz & Stegun 26.7.8 (p.949) */
1959 val = 1./(4.*n);
1960 return norm._cdf(x*(1. - val)/sqrt(1. + x*x*2.*val), lower_tail, log_p);
1961 }
1962
1963 nx = 1 + (x/n)*x;
1964 /* FIXME: This test is probably losing rather than gaining precision,
1965 * now that pbeta(*, log_p = TRUE) is much better.
1966 * Note however that a version of this test *is* needed for x*x > D_MAX */
1967 if(nx > 1e100) { /* <==> x*x > 1e100 * n */
1968 /* Danger of underflow. So use Abramowitz & Stegun 26.5.4
1969 pbeta(z, a, b) ~ z^a(1-z)^b / aB(a,b) ~ z^a / aB(a,b),
1970 with z = 1/nx, a = n/2, b= 1/2 :
1971 */
1972 var lval;
1973 lval = -0.5*n*(2*Math.log(Math.abs(x)) - Math.log(n))
1974 - jstat.logBeta(0.5*n, 0.5) - Math.log(0.5*n);
1975 val = log_p ? lval : Math.exp(lval);
1976 } else {
1977 /*
1978 val = (n > x * x)
1979 // ? pbeta (x * x / (n + x * x), 0.5, n / 2., 0, log_p)
1980 // : pbeta (1. / nx, n / 2., 0.5, 1, log_p);
1981 */
1982 if(n > x * x) {
1983 var beta = new BetaDistribution(0.5, n/2.);
1984 return beta._cdf(x*x/ (n + x * x), false, log_p);
1985 } else {
1986 beta = new BetaDistribution(n / 2., 0.5);
1987 return beta._cdf(1. / nx, true, log_p);
1988 }
1989
1990
1991 }
1992
1993 /* Use "1 - v" if lower_tail and x > 0 (but not both):*/
1994 if(x <= 0.)
1995 lower_tail = !lower_tail;
1996
1997 if(log_p) {
1998 if(lower_tail) return jstat.log1p(-0.5*Math.exp(val));
1999 else return val - M_LN2; /* = log(.5* pbeta(....)) */
2000 }
2001 else {
2002 val /= 2.;
2003 if(lower_tail) {
2004 return (0.5 - val + 0.5);
2005 } else {
2006 return val;
2007 }
2008 }
2009 },
2010 _pnt: function(t, lower_tail, log_p) {
2011
2012 var dof = this._dof;
2013 var ncp = this._mu;
2014 var DT_0, DT_1;
2015
2016 if(lower_tail) {
2017 if(log_p) {
2018 DT_0 = Number.NEGATIVE_INFINITY;
2019 DT_1 = 1.;
2020 } else {
2021 DT_0 = 0.;
2022 DT_1 = 1.;
2023 }
2024 } else {
2025 if(log_p) {
2026 // not lower_tail but log_p
2027 DT_0 = 0.;
2028 DT_1 = Number.NEGATIVE_INFINITY;
2029 } else {
2030 // not lower_tail and not log_p
2031 DT_0 = 1.;
2032 DT_1 = 0.;
2033 }
2034 }
2035
2036 var albeta, a, b, del, errbd, lambda, rxb, tt, x;
2037 var geven, godd, p, q, s, tnc, xeven, xodd;
2038 var it, negdel;
2039
2040 /* note - itrmax and errmax may be changed to suit one's needs. */
2041 var ITRMAX = 1000;
2042 var ERRMAX = 1.e-7;
2043
2044 if(dof <= 0.0) {
2045 return Number.NaN;
2046 } else if (dof == 0.0) {
2047 return this._pt(t);
2048 }
2049
2050 if(!jstat.isFinite(t)) {
2051 return (t < 0) ? DT_0 : DT_1;
2052 }
2053 if(t >= 0.) {
2054 negdel = false;
2055 tt = t;
2056 del = ncp;
2057 } else {
2058 /* We deal quickly with left tail if extreme,
2059 since pt(q, df, ncp) <= pt(0, df, ncp) = \Phi(-ncp) */
2060 if(ncp >= 40 && (!log_p || !lower_tail)) {
2061 return DT_0;
2062 }
2063 negdel = true;
2064 tt = -t;
2065 del = -ncp;
2066 }
2067
2068 if(dof > 4e5 || del*del > 2* Math.LN2 * (-(jstat.DBL_MIN_EXP))) {
2069 /*-- 2nd part: if del > 37.62, then p=0 below
2070 FIXME: test should depend on `df', `tt' AND `del' ! */
2071 /* Approx. from Abramowitz & Stegun 26.7.10 (p.949) */
2072 s=1./(4.*dof);
2073 var norm = new NormalDistribution(del, Math.sqrt(1. + tt*tt*2.*s));
2074 var result = norm._cdf(tt*(1.-s), lower_tail != negdel, log_p);
2075 return result;
2076 }
2077
2078 /* initialize twin series */
2079 /* Guenther, J. (1978). Statist. Computn. Simuln. vol.6, 199. */
2080 x = t * t;
2081 rxb = dof/(x + dof);/* := (1 - x) {x below} -- but more accurately */
2082 x = x / (x + dof);/* in [0,1) */
2083 if (x > 0.) {/* <==> t != 0 */
2084 lambda = del * del;
2085 p = .5 * Math.exp(-.5 * lambda);
2086 if(p == 0.) { // underflow!
2087 console.warn("underflow in _pnt");
2088 return DT_0;
2089 }
2090 q = jstat.SQRT_2dPI * p * del;
2091 s = .5 - p;
2092 if(s < 1e-7) {
2093 s = -0.5 * jstat.expm1(-0.5 * lambda);
2094 }
2095 a = .5;
2096 b = .5 * dof;
2097 /* rxb = (1 - x) ^ b [ ~= 1 - b*x for tiny x --> see 'xeven' below]
2098 * where '(1 - x)' =: rxb {accurately!} above */
2099 rxb = Math.pow(rxb, b);
2100 albeta = jstat.LN_SQRT_PI + jstat.lgamma(b) - jstat.lgamma(.5 + b);
2101 /* TODO: change incompleteBeta function to accept lower_tail and p_log */
2102 xodd = jstat.incompleteBeta(a, b, x);
2103 godd = 2. * rxb * Math.exp(a * Math.log(x) - albeta);
2104 tnc = b * x;
2105 xeven = (tnc < jstat.DBL_EPSILON) ? tnc : 1. - rxb;
2106 geven = tnc * rxb;
2107 tnc = p * xodd + q * xeven;
2108
2109 /* repeat until convergence or iteration limit */
2110 for(it = 1; it <= ITRMAX; it++) {
2111 a += 1.;
2112 xodd -= godd;
2113 xeven -= geven;
2114 godd *= x * (a + b - 1.) / a;
2115 geven *= x * (a + b - .5) / (a + .5);
2116 p *= lambda / (2 * it);
2117 q *= lambda / (2 * it + 1);
2118 tnc += p * xodd + q * xeven;
2119 s -= p;
2120 /* R 2.4.0 added test for rounding error here. */
2121 if(s < -1.e-10) { /* happens e.g. for (t,df,ncp)=(40,10,38.5), after 799 it.*/
2122 //console.write("precision error _pnt");
2123 break;
2124 }
2125 if(s <= 0 && it > 1) break;
2126
2127 errbd = 2. * s * (xodd - godd);
2128
2129 if(Math.abs(errbd) < ERRMAX) break;/*convergence*/
2130 }
2131
2132 if(it == ITRMAX) {
2133 throw "Non-convergence _pnt";
2134 }
2135 } else {
2136 tnc = 0.;
2137 }
2138 norm = new NormalDistribution(0,1);
2139 tnc += norm._cdf(-del, /*lower*/true, /*log_p*/ false);
2140
2141 lower_tail = lower_tail != negdel; /* xor */
2142 if(tnc > 1 - 1e-10 && lower_tail) {
2143 //console.warn("precision error _pnt");
2144 }
2145 var res = jstat.fmin2(tnc, 1.);
2146 if(lower_tail) {
2147 if(log_p) {
2148 return Math.log(res);
2149 } else {
2150 return res;
2151 }
2152 } else {
2153 if(log_p) {
2154 return jstat.log1p(-(res));
2155 } else {
2156 return (0.5 - (res) + 0.5);
2157 }
2158 }
2159 },
2160 getDegreesOfFreedom: function() {
2161 return this._dof;
2162 },
2163 getNonCentralityParameter: function() {
2164 return this._mu;
2165 },
2166 getMean: function() {
2167 if(this._dof > 1) {
2168 var ans = (1/2)*Math.log(this._dof/2) + jstat.lgamma((this._dof-1)/2) - jstat.lgamma(this._dof/2)
2169 return Math.exp(ans) * this._mu;
2170 } else {
2171 return Number.NaN;
2172 }
2173 },
2174 getVariance: function() {
2175 if(this._dof > 2) {
2176 var ans = this._dof * (1 + this._mu*this._mu)/(this._dof-2) - (((this._mu*this._mu * this._dof) / 2) * Math.pow(Math.exp(jstat.lgamma((this._dof - 1)/2)-jstat.lgamma(this._dof/2)), 2));
2177 return ans;
2178 } else {
2179 return Number.NaN;
2180 }
2181 }
2182 });
2183
2184
2185 /******************************************************************************/
2186 /* jQuery Flot graph objects */
2187 /******************************************************************************/
2188
2189 var Plot = Class.extend({
2190 init: function(id, options) {
2191 this._container = '#' + String(id);
2192 this._plots = [];
2193 this._flotObj = null;
2194 this._locked = false;
2195
2196 if(options != null) {
2197 this._options = options;
2198 } else {
2199 this._options = {
2200 };
2201 }
2202
2203 },
2204 getContainer: function() {
2205 return this._container;
2206 },
2207 getGraph: function() {
2208 return this._flotObj;
2209 },
2210 setData: function(data) {
2211 this._plots = data;
2212 },
2213 clear: function() {
2214 this._plots = [];
2215 //this.render();
2216 },
2217 showLegend: function() {
2218 this._options.legend = {
2219 show: true
2220 }
2221 this.render();
2222 },
2223 hideLegend: function() {
2224 this._options.legend = {
2225 show: false
2226 }
2227 this.render();
2228 },
2229 render: function() {
2230 this._flotObj = null;
2231 this._flotObj = $.plot($(this._container), this._plots, this._options);
2232 }
2233 });
2234
2235 var DistributionPlot = Plot.extend({
2236 init: function(id, distribution, range, options) {
2237 this._super(id, options);
2238 this._showPDF = true;
2239 this._showCDF = false;
2240 this._pdfValues = []; // raw values for pdf
2241 this._cdfValues = []; // raw values for cdf
2242 this._maxY = 1;
2243 this._plotType = 'line'; // line, point, both
2244 this._fill = false;
2245 this._distribution = distribution; // underlying PDF
2246 // Range object for the plot
2247 if(range != null && Range.validate(range)) {
2248 this._range = range;
2249 } else {
2250 this._range = this._distribution.getRange(); // no range supplied, use distribution default
2251 }
2252
2253 // render
2254 if(this._distribution != null) {
2255 this._maxY = this._generateValues(); // create the pdf/cdf values in the ctor
2256 } else {
2257 this._options.xaxis = {
2258 min: range.getMinimum(),
2259 max: range.getMaximum()
2260 }
2261 this._options.yaxis = {
2262 max: 1
2263 }
2264 }
2265
2266
2267
2268 this.render();
2269 },
2270 setHover: function(bool) {
2271 if(bool) {
2272 if(this._options.grid == null) {
2273 this._options.grid = {
2274 hoverable: true,
2275 mouseActiveRadius: 25
2276 }
2277 } else {
2278 this._options.grid.hoverable = true,
2279 this._options.grid.mouseActiveRadius = 25
2280 }
2281 function showTooltip(x, y, contents, color) {
2282 $('<div id="jstat_tooltip">' + contents + '</div>').css( {
2283 position: 'absolute',
2284 display: 'none',
2285 top: y + 15,
2286 'font-size': 'small',
2287 left: x + 5,
2288 border: '1px solid ' + color[1],
2289 color: color[2],
2290 padding: '5px',
2291 'background-color': color[0],
2292 opacity: 0.80
2293 }).appendTo("body").show();
2294 }
2295 var previousPoint = null;
2296 $(this._container).bind("plothover", function(event, pos, item) {
2297 $("#x").text(pos.x.toFixed(2));
2298 $("#y").text(pos.y.toFixed(2));
2299 if (item) {
2300 if (previousPoint != item.datapoint) {
2301 previousPoint = item.datapoint;
2302 $("#jstat_tooltip").remove();
2303 var x = jstat.toSigFig(item.datapoint[0],2), y = jstat.toSigFig(item.datapoint[1], 2);
2304 var text = null;
2305 var color = item.series.color;
2306 if(item.series.label == 'PDF') {
2307 text = "P(" + x + ") = " + y;
2308 color = ["#fee", "#fdd", "#C05F5F"];
2309 } else {
2310 // cdf
2311 text = "F(" + x + ") = " + y;
2312 color = ["#eef", "#ddf", "#4A4AC0"];
2313 }
2314 showTooltip(item.pageX, item.pageY, text, color);
2315 }
2316 }
2317 else {
2318 $("#jstat_tooltip").remove();
2319 previousPoint = null;
2320 }
2321 });
2322 $(this._container).bind("mouseleave", function() {
2323 if($('#jstat_tooltip').is(':visible')) {
2324 $('#jstat_tooltip').remove();
2325 previousPoint = null;
2326 }
2327 });
2328 } else {
2329 // unbind
2330 if(this._options.grid == null) {
2331 this._options.grid = {
2332 hoverable: false
2333 }
2334 } else {
2335 this._options.grid.hoverable = false
2336 }
2337 $(this._container).unbind("plothover");
2338 }
2339
2340 this.render();
2341 },
2342 setType: function(type) {
2343 this._plotType = type;
2344 var lines = {};
2345 var points = {};
2346 if(this._plotType == 'line') {
2347 lines.show = true;
2348 points.show = false;
2349 } else if(this._plotType == 'points') {
2350 lines.show = false;
2351 points.show = true;
2352 } else if(this._plotType == 'both') {
2353 lines.show = true;
2354 points.show = true;
2355 }
2356 if(this._options.series == null) {
2357 this._options.series = {
2358 lines: lines,
2359 points: points
2360 }
2361 } else {
2362 if(this._options.series.lines == null) {
2363 this._options.series.lines = lines;
2364 } else {
2365 this._options.series.lines.show = lines.show;
2366 }
2367 if(this._options.series.points == null) {
2368 this._options.series.points = points;
2369 } else {
2370 this._options.series.points.show = points.show;
2371 }
2372 }
2373
2374 this.render();
2375 },
2376 setFill: function(bool) {
2377 this._fill = bool;
2378 if(this._options.series == null) {
2379 this._options.series = {
2380 lines: {
2381 fill: bool
2382 }
2383 }
2384 } else {
2385 if(this._options.series.lines == null) {
2386 this._options.series.lines = {
2387 fill: bool
2388 }
2389 } else {
2390 this._options.series.lines.fill = bool;
2391 }
2392 }
2393 this.render();
2394 },
2395 clear: function() {
2396 this._super();
2397 this._distribution = null;
2398 this._pdfValues = [];
2399 this._cdfValues = [];
2400 this.render();
2401 },
2402 _generateValues: function() {
2403 this._cdfValues = []; // reinitialize the arrays.
2404 this._pdfValues = [];
2405
2406 var xs = this._range.getPoints();
2407
2408 this._options.xaxis = {
2409 min: xs[0],
2410 max: xs[xs.length-1]
2411 }
2412 var pdfs = this._distribution.density(this._range);
2413 var cdfs = this._distribution.cumulativeDensity(this._range);
2414 for(var i = 0; i < xs.length; i++) {
2415 if(pdfs[i] == Number.POSITIVE_INFINITY || pdfs[i] == Number.NEGATIVE_INFINITY) {
2416 pdfs[i] = null;
2417 }
2418 if(cdfs[i] == Number.POSITIVE_INFINITY || cdfs[i] == Number.NEGATIVE_INFINITY) {
2419 cdfs[i] = null;
2420 }
2421 this._pdfValues.push([xs[i], pdfs[i]]);
2422 this._cdfValues.push([xs[i], cdfs[i]]);
2423 }
2424 return jstat.max(pdfs);
2425 },
2426 showPDF: function() {
2427 this._showPDF = true;
2428 this.render();
2429 },
2430 hidePDF: function() {
2431 this._showPDF = false;
2432 this.render();
2433 },
2434 showCDF: function() {
2435 this._showCDF = true;
2436 this.render();
2437 },
2438 hideCDF: function() {
2439 this._showCDF = false;
2440 this.render();
2441 },
2442 setDistribution: function(distribution, range) {
2443 this._distribution = distribution;
2444 if(range != null) {
2445 this._range = range;
2446 } else {
2447 this._range = distribution.getRange();
2448 }
2449 this._maxY = this._generateValues();
2450 this._options.yaxis = {
2451 max: this._maxY*1.1
2452 }
2453
2454 this.render();
2455 },
2456 getDistribution: function() {
2457 return this._distribution;
2458 },
2459 getRange: function() {
2460 return this._range;
2461 },
2462 setRange: function(range) {
2463 this._range = range;
2464 this._generateValues();
2465 this.render();
2466 },
2467 render: function() {
2468 if(this._distribution != null) {
2469 if(this._showPDF && this._showCDF) {
2470 this.setData([{
2471 yaxis: 1,
2472 data:this._pdfValues,
2473 color: 'rgb(237,194,64)',
2474 clickable: false,
2475 hoverable: true,
2476 label: "PDF"
2477 }, {
2478 yaxis: 2,
2479 data:this._cdfValues,
2480 clickable: false,
2481 color: 'rgb(175,216,248)',
2482 hoverable: true,
2483 label: "CDF"
2484 }]);
2485 this._options.yaxis = {
2486 max: this._maxY*1.1
2487 }
2488 } else if(this._showPDF) {
2489 this.setData([{
2490 data:this._pdfValues,
2491 hoverable: true,
2492 color: 'rgb(237,194,64)',
2493 clickable: false,
2494 label: "PDF"
2495 }]);
2496 this._options.yaxis = {
2497 max: this._maxY*1.1
2498 }
2499 } else if(this._showCDF) {
2500 this.setData([{
2501 data:this._cdfValues,
2502 hoverable: true,
2503 color: 'rgb(175,216,248)',
2504 clickable: false,
2505 label: "CDF"
2506 }]);
2507 this._options.yaxis = {
2508 max: 1.1
2509 }
2510 }
2511 } else {
2512 this.setData([]);
2513 }
2514 this._super(); // Call the parent plot method
2515 }
2516 });
2517
2518 var DistributionFactory = {};
2519 DistributionFactory.build = function(json) {
2520 /*
2521 if(json.name == null) {
2522 try{
2523 json = JSON.parse(json);
2524 }
2525 catch(err) {
2526 throw "invalid JSON";
2527 }
2528
2529 // try to parse it
2530 }*/
2531
2532 /*
2533 if(json.name != null) {
2534 var name = json.name;
2535 } else {
2536 throw "Malformed JSON provided to DistributionBuilder " + json;
2537 }
2538 */
2539 if(json.NormalDistribution) {
2540 if(json.NormalDistribution.mean != null && json.NormalDistribution.standardDeviation != null) {
2541 return new NormalDistribution(json.NormalDistribution.mean[0], json.NormalDistribution.standardDeviation[0]);
2542 } else {
2543 throw "Malformed JSON provided to DistributionBuilder " + json;
2544 }
2545 } else if (json.LogNormalDistribution) {
2546 if(json.LogNormalDistribution.location != null && json.LogNormalDistribution.scale != null) {
2547 return new LogNormalDistribution(json.LogNormalDistribution.location[0], json.LogNormalDistribution.scale[0]);
2548 } else {
2549 throw "Malformed JSON provided to DistributionBuilder " + json;
2550 }
2551 } else if (json.BetaDistribution) {
2552 if(json.BetaDistribution.alpha != null && json.BetaDistribution.beta != null) {
2553 return new BetaDistribution(json.BetaDistribution.alpha[0], json.BetaDistribution.beta[0]);
2554 } else {
2555 throw "Malformed JSON provided to DistributionBuilder " + json;
2556 }
2557 } else if (json.GammaDistribution) {
2558 if(json.GammaDistribution.shape != null && json.GammaDistribution.scale != null) {
2559 return new GammaDistribution(json.GammaDistribution.shape[0], json.GammaDistribution.scale[0]);
2560 } else {
2561 throw "Malformed JSON provided to DistributionBuilder " + json;
2562 }
2563 } else if (json.StudentTDistribution) {
2564 if(json.StudentTDistribution.degreesOfFreedom != null && json.StudentTDistribution.nonCentralityParameter != null) {
2565 return new StudentTDistribution(json.StudentTDistribution.degreesOfFreedom[0], json.StudentTDistribution.nonCentralityParameter[0]);
2566 } else if(json.StudentTDistribution.degreesOfFreedom != null) {
2567 return new StudentTDistribution(json.StudentTDistribution.degreesOfFreedom[0]);
2568 } else {
2569 throw "Malformed JSON provided to DistributionBuilder " + json;
2570 }
2571 } else {
2572 throw "Malformed JSON provided to DistributionBuilder " + json;
2573 }
2574 }
2575
2576