Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
JavaScript library for calculating critical values and upper probabilities of common statistical distributions
/*
* NAME
*
* statistics-distributions.js - JavaScript library for calculating
* critical values and upper probabilities of common statistical
* distributions
*
* SYNOPSIS
*
*
* // Chi-squared-crit (2 degrees of freedom, 95th percentile = 0.05 level
* chisqrdistr(2, .05)
*
* // u-crit (95th percentile = 0.05 level)
* udistr(.05);
*
* // t-crit (1 degree of freedom, 99.5th percentile = 0.005 level)
* tdistr(1,.005);
*
* // F-crit (1 degree of freedom in numerator, 3 degrees of freedom
* // in denominator, 99th percentile = 0.01 level)
* fdistr(1,3,.01);
*
* // upper probability of the u distribution (u = -0.85): Q(u) = 1-G(u)
* uprob(-0.85);
*
* // upper probability of the chi-square distribution
* // (3 degrees of freedom, chi-squared = 6.25): Q = 1-G
* chisqrprob(3,6.25);
*
* // upper probability of the t distribution
* // (3 degrees of freedom, t = 6.251): Q = 1-G
* tprob(3,6.251);
*
* // upper probability of the F distribution
* // (3 degrees of freedom in numerator, 5 degrees of freedom in
* // denominator, F = 6.25): Q = 1-G
* fprob(3,5,.625);
*
*
* DESCRIPTION
*
* This library calculates percentage points (5 significant digits) of the u
* (standard normal) distribution, the student's t distribution, the
* chi-square distribution and the F distribution. It can also calculate the
* upper probability (5 significant digits) of the u (standard normal), the
* chi-square, the t and the F distribution.
*
* These critical values are needed to perform statistical tests, like the u
* test, the t test, the F test and the chi-squared test, and to calculate
* confidence intervals.
*
* If you are interested in more precise algorithms you could look at:
* StatLib: http://lib.stat.cmu.edu/apstat/ ;
* Applied Statistics Algorithms by Griffiths, P. and Hill, I.D.
* , Ellis Horwood: Chichester (1985)
*
* BUGS
*
* This port was produced from the Perl module Statistics::Distributions
* that has had no bug reports in several years. If you find a bug then
* please double-check that JavaScript does not thing the numbers you are
* passing in are strings. (You can subtract 0 from them as you pass them
* in so that "5" is properly understood to be 5.) If you have passed in a
* number then please contact the author
*
* AUTHOR
*
* Ben Tilly <btilly@gmail.com>
*
* Originl Perl version by Michael Kospach <mike.perl@gmx.at>
*
* Nice formating, simplification and bug repair by Matthias Trautner Kromann
* <mtk@id.cbs.dk>
*
* COPYRIGHT
*
* Copyright 2008 Ben Tilly.
*
* This library is free software; you can redistribute it and/or modify it
* under the same terms as Perl itself. This means under either the Perl
* Artistic License or the GPL v1 or later.
*/
var SIGNIFICANT = 5; // number of significant digits to be returned
function chisqrdistr ($n, $p) {
if ($n <= 0 || Math.abs($n) - Math.abs(integer($n)) != 0) {
throw("Invalid n: $n\n"); /* degree of freedom */
}
if ($p <= 0 || $p > 1) {
throw("Invalid p: $p\n");
}
return precision_string(_subchisqr($n-0, $p-0));
}
function udistr ($p) {
if ($p > 1 || $p <= 0) {
throw("Invalid p: $p\n");
}
return precision_string(_subu($p-0));
}
function tdistr ($n, $p) {
if ($n <= 0 || Math.abs($n) - Math.abs(integer($n)) != 0) {
throw("Invalid n: $n\n");
}
if ($p <= 0 || $p >= 1) {
throw("Invalid p: $p\n");
}
return precision_string(_subt($n-0, $p-0));
}
function fdistr ($n, $m, $p) {
if (($n<=0) || ((Math.abs($n)-(Math.abs(integer($n))))!=0)) {
throw("Invalid n: $n\n"); /* first degree of freedom */
}
if (($m<=0) || ((Math.abs($m)-(Math.abs(integer($m))))!=0)) {
throw("Invalid m: $m\n"); /* second degree of freedom */
}
if (($p<=0) || ($p>1)) {
throw("Invalid p: $p\n");
}
return precision_string(_subf($n-0, $m-0, $p-0));
}
function uprob ($x) {
return precision_string(_subuprob($x-0));
}
function chisqrprob ($n,$x) {
if (($n <= 0) || ((Math.abs($n) - (Math.abs(integer($n)))) != 0)) {
throw("Invalid n: $n\n"); /* degree of freedom */
}
return precision_string(_subchisqrprob($n-0, $x-0));
}
function tprob ($n, $x) {
if (($n <= 0) || ((Math.abs($n) - Math.abs(integer($n))) !=0)) {
throw("Invalid n: $n\n"); /* degree of freedom */
}
return precision_string(_subtprob($n-0, $x-0));
}
function fprob ($n, $m, $x) {
if (($n<=0) || ((Math.abs($n)-(Math.abs(integer($n))))!=0)) {
throw("Invalid n: $n\n"); /* first degree of freedom */
}
if (($m<=0) || ((Math.abs($m)-(Math.abs(integer($m))))!=0)) {
throw("Invalid m: $m\n"); /* second degree of freedom */
}
return precision_string(_subfprob($n-0, $m-0, $x-0));
}
function _subfprob ($n, $m, $x) {
var $p;
if ($x<=0) {
$p=1;
} else if ($m % 2 == 0) {
var $z = $m / ($m + $n * $x);
var $a = 1;
for (var $i = $m - 2; $i >= 2; $i -= 2) {
$a = 1 + ($n + $i - 2) / $i * $z * $a;
}
$p = 1 - Math.pow((1 - $z), ($n / 2) * $a);
} else if ($n % 2 == 0) {
var $z = $n * $x / ($m + $n * $x);
var $a = 1;
for (var $i = $n - 2; $i >= 2; $i -= 2) {
$a = 1 + ($m + $i - 2) / $i * $z * $a;
}
$p = Math.pow((1 - $z), ($m / 2)) * $a;
} else {
var $y = Math.atan2(Math.sqrt($n * $x / $m), 1);
var $z = Math.pow(Math.sin($y), 2);
var $a = ($n == 1) ? 0 : 1;
for (var $i = $n - 2; $i >= 3; $i -= 2) {
$a = 1 + ($m + $i - 2) / $i * $z * $a;
}
var $b = Math.PI;
for (var $i = 2; $i <= $m - 1; $i += 2) {
$b *= ($i - 1) / $i;
}
var $p1 = 2 / $b * Math.sin($y) * Math.pow(Math.cos($y), $m) * $a;
$z = Math.pow(Math.cos($y), 2);
$a = ($m == 1) ? 0 : 1;
for (var $i = $m-2; $i >= 3; $i -= 2) {
$a = 1 + ($i - 1) / $i * $z * $a;
}
$p = max(0, $p1 + 1 - 2 * $y / Math.PI
- 2 / Math.PI * Math.sin($y) * Math.cos($y) * $a);
}
return $p;
}
function _subchisqrprob ($n,$x) {
var $p;
if ($x <= 0) {
$p = 1;
} else if ($n > 100) {
$p = _subuprob((Math.pow(($x / $n), 1/3)
- (1 - 2/9/$n)) / Math.sqrt(2/9/$n));
} else if ($x > 400) {
$p = 0;
} else {
var $a;
var $i;
var $i1;
if (($n % 2) != 0) {
$p = 2 * _subuprob(Math.sqrt($x));
$a = Math.sqrt(2/Math.PI) * Math.exp(-$x/2) / Math.sqrt($x);
$i1 = 1;
} else {
$p = $a = Math.exp(-$x/2);
$i1 = 2;
}
for ($i = $i1; $i <= ($n-2); $i += 2) {
$a *= $x / $i;
$p += $a;
}
}
return $p;
}
function _subu ($p) {
var $y = -Math.log(4 * $p * (1 - $p));
var $x = Math.sqrt(
$y * (1.570796288
+ $y * (.03706987906
+ $y * (-.8364353589E-3
+ $y *(-.2250947176E-3
+ $y * (.6841218299E-5
+ $y * (0.5824238515E-5
+ $y * (-.104527497E-5
+ $y * (.8360937017E-7
+ $y * (-.3231081277E-8
+ $y * (.3657763036E-10
+ $y *.6936233982E-12)))))))))));
if ($p>.5)
$x = -$x;
return $x;
}
function _subuprob ($x) {
var $p = 0; /* if ($absx > 100) */
var $absx = Math.abs($x);
if ($absx < 1.9) {
$p = Math.pow((1 +
$absx * (.049867347
+ $absx * (.0211410061
+ $absx * (.0032776263
+ $absx * (.0000380036
+ $absx * (.0000488906
+ $absx * .000005383)))))), -16)/2;
} else if ($absx <= 100) {
for (var $i = 18; $i >= 1; $i--) {
$p = $i / ($absx + $p);
}
$p = Math.exp(-.5 * $absx * $absx)
/ Math.sqrt(2 * Math.PI) / ($absx + $p);
}
if ($x<0)
$p = 1 - $p;
return $p;
}
function _subt ($n, $p) {
if ($p >= 1 || $p <= 0) {
throw("Invalid p: $p\n");
}
if ($p == 0.5) {
return 0;
} else if ($p < 0.5) {
return - _subt($n, 1 - $p);
}
var $u = _subu($p);
var $u2 = Math.pow($u, 2);
var $a = ($u2 + 1) / 4;
var $b = ((5 * $u2 + 16) * $u2 + 3) / 96;
var $c = (((3 * $u2 + 19) * $u2 + 17) * $u2 - 15) / 384;
var $d = ((((79 * $u2 + 776) * $u2 + 1482) * $u2 - 1920) * $u2 - 945)
/ 92160;
var $e = (((((27 * $u2 + 339) * $u2 + 930) * $u2 - 1782) * $u2 - 765) * $u2
+ 17955) / 368640;
var $x = $u * (1 + ($a + ($b + ($c + ($d + $e / $n) / $n) / $n) / $n) / $n);
if ($n <= Math.pow(log10($p), 2) + 3) {
var $round;
do {
var $p1 = _subtprob($n, $x);
var $n1 = $n + 1;
var $delta = ($p1 - $p)
/ Math.exp(($n1 * Math.log($n1 / ($n + $x * $x))
+ Math.log($n/$n1/2/Math.PI) - 1
+ (1/$n1 - 1/$n) / 6) / 2);
$x += $delta;
$round = round_to_precision($delta, Math.abs(integer(log10(Math.abs($x))-4)));
} while (($x) && ($round != 0));
}
return $x;
}
function _subtprob ($n, $x) {
var $a;
var $b;
var $w = Math.atan2($x / Math.sqrt($n), 1);
var $z = Math.pow(Math.cos($w), 2);
var $y = 1;
for (var $i = $n-2; $i >= 2; $i -= 2) {
$y = 1 + ($i-1) / $i * $z * $y;
}
if ($n % 2 == 0) {
$a = Math.sin($w)/2;
$b = .5;
} else {
$a = ($n == 1) ? 0 : Math.sin($w)*Math.cos($w)/Math.PI;
$b= .5 + $w/Math.PI;
}
return max(0, 1 - $b - $a * $y);
}
function _subf ($n, $m, $p) {
var $x;
if ($p >= 1 || $p <= 0) {
throw("Invalid p: $p\n");
}
if ($p == 1) {
$x = 0;
} else if ($m == 1) {
$x = 1 / Math.pow(_subt($n, 0.5 - $p / 2), 2);
} else if ($n == 1) {
$x = Math.pow(_subt($m, $p/2), 2);
} else if ($m == 2) {
var $u = _subchisqr($m, 1 - $p);
var $a = $m - 2;
$x = 1 / ($u / $m * (1 +
(($u - $a) / 2 +
(((4 * $u - 11 * $a) * $u + $a * (7 * $m - 10)) / 24 +
(((2 * $u - 10 * $a) * $u + $a * (17 * $m - 26)) * $u
- $a * $a * (9 * $m - 6)
)/48/$n
)/$n
)/$n));
} else if ($n > $m) {
$x = 1 / _subf2($m, $n, 1 - $p)
} else {
$x = _subf2($n, $m, $p)
}
return $x;
}
function _subf2 ($n, $m, $p) {
var $u = _subchisqr($n, $p);
var $n2 = $n - 2;
var $x = $u / $n *
(1 +
(($u - $n2) / 2 +
(((4 * $u - 11 * $n2) * $u + $n2 * (7 * $n - 10)) / 24 +
(((2 * $u - 10 * $n2) * $u + $n2 * (17 * $n - 26)) * $u
- $n2 * $n2 * (9 * $n - 6)) / 48 / $m) / $m) / $m);
var $delta;
do {
var $z = Math.exp(
(($n+$m) * Math.log(($n+$m) / ($n * $x + $m))
+ ($n - 2) * Math.log($x)
+ Math.log($n * $m / ($n+$m))
- Math.log(4 * Math.PI)
- (1/$n + 1/$m - 1/($n+$m))/6
)/2);
$delta = (_subfprob($n, $m, $x) - $p) / $z;
$x += $delta;
} while (Math.abs($delta)>3e-4);
return $x;
}
function _subchisqr ($n, $p) {
var $x;
if (($p > 1) || ($p <= 0)) {
throw("Invalid p: $p\n");
} else if ($p == 1){
$x = 0;
} else if ($n == 1) {
$x = Math.pow(_subu($p / 2), 2);
} else if ($n == 2) {
$x = -2 * Math.log($p);
} else {
var $u = _subu($p);
var $u2 = $u * $u;
$x = max(0, $n + Math.sqrt(2 * $n) * $u
+ 2/3 * ($u2 - 1)
+ $u * ($u2 - 7) / 9 / Math.sqrt(2 * $n)
- 2/405 / $n * ($u2 * (3 *$u2 + 7) - 16));
if ($n <= 100) {
var $x0;
var $p1;
var $z;
do {
$x0 = $x;
if ($x < 0) {
$p1 = 1;
} else if ($n>100) {
$p1 = _subuprob((Math.pow(($x / $n), (1/3)) - (1 - 2/9/$n))
/ Math.sqrt(2/9/$n));
} else if ($x>400) {
$p1 = 0;
} else {
var $i0
var $a;
if (($n % 2) != 0) {
$p1 = 2 * _subuprob(Math.sqrt($x));
$a = Math.sqrt(2/Math.PI) * Math.exp(-$x/2) / Math.sqrt($x);
$i0 = 1;
} else {
$p1 = $a = Math.exp(-$x/2);
$i0 = 2;
}
for (var $i = $i0; $i <= $n-2; $i += 2) {
$a *= $x / $i;
$p1 += $a;
}
}
$z = Math.exp((($n-1) * Math.log($x/$n) - Math.log(4*Math.PI*$x)
+ $n - $x - 1/$n/6) / 2);
$x += ($p1 - $p) / $z;
$x = round_to_precision($x, 5);
} while (($n < 31) && (Math.abs($x0 - $x) > 1e-4));
}
}
return $x;
}
function log10 ($n) {
return Math.log($n) / Math.log(10);
}
function max () {
var $max = arguments[0];
for (var $i = 0; i < arguments.length; i++) {
if ($max < arguments[$i])
$max = arguments[$i];
}
return $max;
}
function min () {
var $min = arguments[0];
for (var $i = 0; i < arguments.length; i++) {
if ($min > arguments[$i])
$min = arguments[$i];
}
return $min;
}
function precision ($x) {
return Math.abs(integer(log10(Math.abs($x)) - SIGNIFICANT));
}
function precision_string ($x) {
if ($x) {
return round_to_precision($x, precision($x));
} else {
return "0";
}
}
function round_to_precision ($x, $p) {
$x = $x * Math.pow(10, $p);
$x = Math.round($x);
return $x / Math.pow(10, $p);
}
function integer ($i) {
if ($i > 0)
return Math.floor($i);
else
return Math.ceil($i);
}
@slfan2013

This comment has been minimized.

Copy link

@slfan2013 slfan2013 commented Aug 20, 2016

Hi
Why the tprob(5,0) returns 0?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment