|
0
|
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
|
