[system] / trunk / pg / lib / Distributions.pm Repository:
ViewVC logotype

View of /trunk/pg/lib/Distributions.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6346 - (download) (as text) (annotate)
Sat Jul 10 12:39:40 2010 UTC (9 years, 7 months ago) by gage
File size: 11758 byte(s)
Merging changes gage branch  gage_dev/pg

removed dependence on AUTOLOAD	which does not work well with newer versions of Safe.pm.  It wasn't needed 
in any case.  There remain other incompatibilies of WeBWorK with Safe.pm 2.27

Added more support for WARN_MESSAGE  and DEBUG_MESSAGE

Changed List.pm to ChoiceList.pm  to remove confusion with MathObjects List object

Additional support for geogebra applets



    1 #package Statistics::Distributions;
    2 package Distributions;
    3 
    4 # In order to use the functions from this module in your WebWork problems, you have to load PGstatisticsmacros.pl
    5 
    6 # The documentation at the end of this file is slightly modified for WebWork by Maria Voloshina, University of Rochester.
    7 
    8 use strict;
    9 use vars qw($VERSION @ISA @EXPORT @EXPORT_OK);
   10 use constant PI => 3.1415926536;
   11 use constant SIGNIFICANT => 6; # number of significant digits to be returned
   12 
   13 require Exporter;
   14 
   15 #@ISA = qw(Exporter AutoLoader);
   16 @ISA = qw(Exporter);
   17 # Items to export into callers namespace by default. Note: do not export
   18 # names by default without a very good reason. Use EXPORT_OK instead.
   19 # Do not simply export all your public functions/methods/constants.
   20 @EXPORT_OK = qw(chisqrdistr tdistr fdistr udistr uprob chisqrprob tprob fprob);
   21 $VERSION = '0.07';
   22 
   23 # Preloaded methods go here.
   24 
   25 sub chisqrdistr { # Percentage points  X^2(x^2,n)
   26   my ($n, $p) = @_;
   27   if ($n <= 0 || abs($n) - abs(int($n)) != 0) {
   28     die "Invalid n: $n\n"; # degree of freedom
   29   }
   30   if ($p <= 0 || $p > 1) {
   31     die "Invalid p: $p\n";
   32   }
   33   return precision_string(_subchisqr($n, $p));
   34 }
   35 
   36 sub udistr { # Percentage points   N(0,1^2)
   37   my ($p) = (@_);
   38   if ($p > 1 || $p <= 0) {
   39     die "Invalid p: $p\n";
   40   }
   41   return precision_string(_subu($p));
   42 }
   43 
   44 sub tdistr { # Percentage points   t(x,n)
   45   my ($n, $p) = @_;
   46   if ($n <= 0 || abs($n) - abs(int($n)) != 0) {
   47     die "Invalid n: $n\n";
   48   }
   49   if ($p <= 0 || $p >= 1) {
   50     die "Invalid p: $p\n";
   51   }
   52   return precision_string(_subt($n, $p));
   53 }
   54 
   55 sub fdistr { # Percentage points  F(x,n1,n2)
   56   my ($n, $m, $p) = @_;
   57   if (($n<=0) || ((abs($n)-(abs(int($n))))!=0)) {
   58     die "Invalid n: $n\n"; # first degree of freedom
   59   }
   60   if (($m<=0) || ((abs($m)-(abs(int($m))))!=0)) {
   61     die "Invalid m: $m\n"; # second degree of freedom
   62   }
   63   if (($p<=0) || ($p>1)) {
   64     die "Invalid p: $p\n";
   65   }
   66   return precision_string(_subf($n, $m, $p));
   67 }
   68 
   69 sub uprob { # Upper probability   N(0,1^2)
   70   my ($x) = @_;
   71   return precision_string(_subuprob($x));
   72 }
   73 
   74 sub chisqrprob { # Upper probability   X^2(x^2,n)
   75   my ($n,$x) = @_;
   76   if (($n <= 0) || ((abs($n) - (abs(int($n)))) != 0)) {
   77     die "Invalid n: $n\n"; # degree of freedom
   78   }
   79   return precision_string(_subchisqrprob($n, $x));
   80 }
   81 
   82 sub tprob { # Upper probability   t(x,n)
   83   my ($n, $x) = @_;
   84   if (($n <= 0) || ((abs($n) - abs(int($n))) !=0)) {
   85     die "Invalid n: $n\n"; # degree of freedom
   86   }
   87   return precision_string(_subtprob($n, $x));
   88 }
   89 
   90 sub fprob { # Upper probability   F(x,n1,n2)
   91   my ($n, $m, $x) = @_;
   92   if (($n<=0) || ((abs($n)-(abs(int($n))))!=0)) {
   93     die "Invalid n: $n\n"; # first degree of freedom
   94   }
   95   if (($m<=0) || ((abs($m)-(abs(int($m))))!=0)) {
   96     die "Invalid m: $m\n"; # second degree of freedom
   97   }
   98   return precision_string(_subfprob($n, $m, $x));
   99 }
  100 
  101 
  102 sub _subfprob {
  103   my ($n, $m, $x) = @_;
  104   my $p;
  105 
  106   if ($x<=0) {
  107     $p=1;
  108   } elsif ($m % 2 == 0) {
  109     my $z = $m / ($m + $n * $x);
  110     my $a = 1;
  111     for (my $i = $m - 2; $i >= 2; $i -= 2) {
  112       $a = 1 + ($n + $i - 2) / $i * $z * $a;
  113     }
  114     $p = 1 - ((1 - $z) ** ($n / 2) * $a);
  115   } elsif ($n % 2 == 0) {
  116     my $z = $n * $x / ($m + $n * $x);
  117     my $a = 1;
  118     for (my $i = $n - 2; $i >= 2; $i -= 2) {
  119       $a = 1 + ($m + $i - 2) / $i * $z * $a;
  120     }
  121     $p = (1 - $z) ** ($m / 2) * $a;
  122   } else {
  123     my $y = atan2(sqrt($n * $x / $m), 1);
  124     my $z = sin($y) ** 2;
  125     my $a = ($n == 1) ? 0 : 1;
  126     for (my $i = $n - 2; $i >= 3; $i -= 2) {
  127       $a = 1 + ($m + $i - 2) / $i * $z * $a;
  128     }
  129     my $b = PI;
  130     for (my $i = 2; $i <= $m - 1; $i += 2) {
  131       $b *= ($i - 1) / $i;
  132     }
  133     my $p1 = 2 / $b * sin($y) * cos($y) ** $m * $a;
  134 
  135     $z = cos($y) ** 2;
  136     $a = ($m == 1) ? 0 : 1;
  137     for (my $i = $m-2; $i >= 3; $i -= 2) {
  138       $a = 1 + ($i - 1) / $i * $z * $a;
  139     }
  140     $p = max(0, $p1 + 1 - 2 * $y / PI
  141       - 2 / PI * sin($y) * cos($y) * $a);
  142   }
  143   return $p;
  144 }
  145 
  146 
  147 sub _subchisqrprob {
  148   my ($n,$x) = @_;
  149   my $p;
  150 
  151   if ($x <= 0) {
  152     $p = 1;
  153   } elsif ($n > 100) {
  154     $p = _subuprob((($x / $n) ** (1/3)
  155         - (1 - 2/9/$n)) / sqrt(2/9/$n));
  156   } elsif ($x > 400) {
  157     $p = 0;
  158   } else {
  159     my ($a, $i, $i1);
  160     if (($n % 2) != 0) {
  161       $p = 2 * _subuprob(sqrt($x));
  162       $a = sqrt(2/PI) * exp(-$x/2) / sqrt($x);
  163       $i1 = 1;
  164     } else {
  165       $p = $a = exp(-$x/2);
  166       $i1 = 2;
  167     }
  168 
  169     for ($i = $i1; $i <= ($n-2); $i += 2) {
  170       $a *= $x / $i;
  171       $p += $a;
  172     }
  173   }
  174   return $p;
  175 }
  176 
  177 sub _subu {
  178   my ($p) = @_;
  179   my $y = -log(4 * $p * (1 - $p));
  180   my $x = sqrt(
  181     $y * (1.570796288
  182       + $y * (.03706987906
  183         + $y * (-.8364353589E-3
  184         + $y *(-.2250947176E-3
  185           + $y * (.6841218299E-5
  186           + $y * (0.5824238515E-5
  187           + $y * (-.104527497E-5
  188             + $y * (.8360937017E-7
  189             + $y * (-.3231081277E-8
  190               + $y * (.3657763036E-10
  191               + $y *.6936233982E-12)))))))))));
  192   $x = -$x if ($p>.5);
  193   return $x;
  194 }
  195 
  196 sub _subuprob {
  197   my ($x) = @_;
  198   my $p = 0; # if ($absx > 100)
  199   my $absx = abs($x);
  200 
  201   if ($absx < 1.9) {
  202     $p = (1 +
  203       $absx * (.049867347
  204         + $absx * (.0211410061
  205           + $absx * (.0032776263
  206           + $absx * (.0000380036
  207           + $absx * (.0000488906
  208             + $absx * .000005383)))))) ** -16/2;
  209   } elsif ($absx <= 100) {
  210     for (my $i = 18; $i >= 1; $i--) {
  211       $p = $i / ($absx + $p);
  212     }
  213     $p = exp(-.5 * $absx * $absx)
  214       / sqrt(2 * PI) / ($absx + $p);
  215   }
  216 
  217   $p = 1 - $p if ($x<0);
  218   return $p;
  219 }
  220 
  221 
  222 sub _subt {
  223   my ($n, $p) = @_;
  224 
  225   if ($p >= 1 || $p <= 0) {
  226     die "Invalid p: $p\n";
  227   }
  228 
  229   if ($p == 0.5) {
  230     return 0;
  231   } elsif ($p < 0.5) {
  232     return - _subt($n, 1 - $p);
  233   }
  234 
  235   my $u = _subu($p);
  236   my $u2 = $u ** 2;
  237 
  238   my $a = ($u2 + 1) / 4;
  239   my $b = ((5 * $u2 + 16) * $u2 + 3) / 96;
  240   my $c = (((3 * $u2 + 19) * $u2 + 17) * $u2 - 15) / 384;
  241   my $d = ((((79 * $u2 + 776) * $u2 + 1482) * $u2 - 1920) * $u2 - 945)
  242         / 92160;
  243   my $e = (((((27 * $u2 + 339) * $u2 + 930) * $u2 - 1782) * $u2 - 765) * $u2
  244       + 17955) / 368640;
  245 
  246   my $x = $u * (1 + ($a + ($b + ($c + ($d + $e / $n) / $n) / $n) / $n) / $n);
  247 
  248   if ($n <= log10($p) ** 2 + 3) {
  249     my $round;
  250     do {
  251       my $p1 = _subtprob($n, $x);
  252       my $n1 = $n + 1;
  253       my $delta = ($p1 - $p)
  254         / exp(($n1 * log($n1 / ($n + $x * $x))
  255           + log($n/$n1/2/PI) - 1
  256           + (1/$n1 - 1/$n) / 6) / 2);
  257       $x += $delta;
  258       $round = sprintf("%.".abs(int(log10(abs $x)-4))."f",$delta);
  259     } while (($x) && ($round != 0));
  260   }
  261   return $x;
  262 }
  263 
  264 sub _subtprob {
  265   my ($n, $x) = @_;
  266 
  267   my ($a,$b);
  268   my $w = atan2($x / sqrt($n), 1);
  269   my $z = cos($w) ** 2;
  270   my $y = 1;
  271 
  272   for (my $i = $n-2; $i >= 2; $i -= 2) {
  273     $y = 1 + ($i-1) / $i * $z * $y;
  274   }
  275 
  276   if ($n % 2 == 0) {
  277     $a = sin($w)/2;
  278     $b = .5;
  279   } else {
  280     $a = ($n == 1) ? 0 : sin($w)*cos($w)/PI;
  281     $b= .5 + $w/PI;
  282   }
  283   return max(0, 1 - $b - $a * $y);
  284 }
  285 
  286 sub _subf {
  287   my ($n, $m, $p) = @_;
  288   my $x;
  289 
  290   if ($p >= 1 || $p <= 0) {
  291     die "Invalid p: $p\n";
  292   }
  293 
  294   if ($p == 1) {
  295     $x = 0;
  296   } elsif ($m == 1) {
  297     $x = 1 / (_subt($n, 0.5 - $p / 2) ** 2);
  298   } elsif ($n == 1) {
  299     $x = _subt($m, $p/2) ** 2;
  300   } elsif ($m == 2) {
  301     my $u = _subchisqr($m, 1 - $p);
  302     my $a = $m - 2;
  303     $x = 1 / ($u / $m * (1 +
  304       (($u - $a) / 2 +
  305         (((4 * $u - 11 * $a) * $u + $a * (7 * $m - 10)) / 24 +
  306           (((2 * $u - 10 * $a) * $u + $a * (17 * $m - 26)) * $u
  307             - $a * $a * (9 * $m - 6)
  308           )/48/$n
  309         )/$n
  310       )/$n));
  311   } elsif ($n > $m) {
  312     $x = 1 / _subf2($m, $n, 1 - $p)
  313   } else {
  314     $x = _subf2($n, $m, $p)
  315   }
  316   return $x;
  317 }
  318 
  319 sub _subf2 {
  320   my ($n, $m, $p) = @_;
  321   my $u = _subchisqr($n, $p);
  322   my $n2 = $n - 2;
  323   my $x = $u / $n *
  324     (1 +
  325       (($u - $n2) / 2 +
  326         (((4 * $u - 11 * $n2) * $u + $n2 * (7 * $n - 10)) / 24 +
  327           (((2 * $u - 10 * $n2) * $u + $n2 * (17 * $n - 26)) * $u
  328             - $n2 * $n2 * (9 * $n - 6)) / 48 / $m) / $m) / $m);
  329   my $delta;
  330   do {
  331     my $z = exp(
  332       (($n+$m) * log(($n+$m) / ($n * $x + $m))
  333         + ($n - 2) * log($x)
  334         + log($n * $m / ($n+$m))
  335         - log(4 * PI)
  336         - (1/$n  + 1/$m - 1/($n+$m))/6
  337       )/2);
  338     $delta = (_subfprob($n, $m, $x) - $p) / $z;
  339     $x += $delta;
  340   } while (abs($delta)>3e-4);
  341   return $x;
  342 }
  343 
  344 sub _subchisqr {
  345   my ($n, $p) = @_;
  346   my $x;
  347 
  348   if (($p > 1) || ($p <= 0)) {
  349     die "Invalid p: $p\n";
  350   } elsif ($p == 1){
  351     $x = 0;
  352   } elsif ($n == 1) {
  353     $x = _subu($p / 2) ** 2;
  354   } elsif ($n == 2) {
  355     $x = -2 * log($p);
  356   } else {
  357     my $u = _subu($p);
  358     my $u2 = $u * $u;
  359 
  360     $x = max(0, $n + sqrt(2 * $n) * $u
  361       + 2/3 * ($u2 - 1)
  362       + $u * ($u2 - 7) / 9 / sqrt(2 * $n)
  363       - 2/405 / $n * ($u2 * (3 *$u2 + 7) - 16));
  364 
  365     if ($n <= 100) {
  366       my ($x0, $p1, $z);
  367       do {
  368         $x0 = $x;
  369         if ($x < 0) {
  370           $p1 = 1;
  371         } elsif ($n>100) {
  372           $p1 = _subuprob((($x / $n)**(1/3) - (1 - 2/9/$n))
  373             / sqrt(2/9/$n));
  374         } elsif ($x>400) {
  375           $p1 = 0;
  376         } else {
  377           my ($i0, $a);
  378           if (($n % 2) != 0) {
  379             $p1 = 2 * _subuprob(sqrt($x));
  380             $a = sqrt(2/PI) * exp(-$x/2) / sqrt($x);
  381             $i0 = 1;
  382           } else {
  383             $p1 = $a = exp(-$x/2);
  384             $i0 = 2;
  385           }
  386 
  387           for (my $i = $i0; $i <= $n-2; $i += 2) {
  388             $a *= $x / $i;
  389             $p1 += $a;
  390           }
  391         }
  392         $z = exp((($n-1) * log($x/$n) - log(4*PI*$x)
  393           + $n - $x - 1/$n/6) / 2);
  394         $x += ($p1 - $p) / $z;
  395         $x = sprintf("%.5f", $x);
  396       } while (($n < 31) && (abs($x0 - $x) > 1e-4));
  397     }
  398   }
  399   return $x;
  400 }
  401 
  402 sub log10 {
  403   my $n = shift;
  404   return log($n) / log(10);
  405 }
  406 
  407 sub max {
  408   my $max = shift;
  409   my $next;
  410   while (@_) {
  411     $next = shift;
  412     $max = $next if ($next > $max);
  413   }
  414   return $max;
  415 }
  416 
  417 sub min {
  418   my $min = shift;
  419   my $next;
  420   while (@_) {
  421     $next = shift;
  422     $min = $next if ($next < $min);
  423   }
  424   return $min;
  425 }
  426 
  427 sub precision {
  428   my ($x) = @_;
  429   return abs int(log10(abs $x) - SIGNIFICANT);
  430 }
  431 
  432 sub precision_string {
  433   my ($x) = @_;
  434   if ($x) {
  435     return sprintf "%." . precision($x) . "f", $x;
  436   } else {
  437     return "0";
  438   }
  439 }
  440 
  441 
  442 # Autoload methods go after =cut, and are processed by the autosplit program.
  443 
  444 1;
  445 __END__
  446 
  447 # Below is the stub of documentation for your module.
  448 
  449 =head1 NAME
  450 
  451 Statistics::Distributions - Perl module for calculating probabilities and critical values of common statistical distributions
  452 
  453 =head1 SYNOPSIS
  454 
  455   use Statistics::Distributions;
  456 
  457   $uprob=uprob(-0.85);
  458   print "upper probability of the u distribution: Q(u) = 1-G(u) (u=1.43) = $uprob";
  459 
  460   $tprob=tprob(3, 6.251);
  461   print "upper probability of the t distribution: Q = 1-G (3 degrees of freedom  , t = 6.251) = $tprob";
  462 
  463   $chisprob=chisqrprob(3, 6.25);
  464   print "upper probability of the chi-square distribution: Q = 1-G (3 degrees of freedom, chi-squared = 6.25) = $chisprob";
  465 
  466   $fprob=fprob(3,5, .625);
  467   print "upper probability of the F distribution: Q = 1-G (3 degrees of freedom  in numerator, 5 degrees of freedom in denominator, F $
  468 
  469   $u=udistr(.05);
  470   print "u-crit (95th percentile = 0.05 level) = $u";
  471 
  472   $t=tdistr(1, .005);
  473   print "t-crit (1 degree of freedom, 99.5th percentile = 0.005 level) =$t";
  474 
  475   $chis=chisqrdistr(2, .05);
  476   print "Chi-squared-crit (2 degrees of freedom, 95th percentile = 0.05 level) = $chis";
  477 
  478   $f=fdistr(1,3, .01);
  479   print "F-crit (1 degree of freedom in numerator, 3 degrees of freedom in denominator, 99th percentile = 0.01 level) = $f";
  480 
  481 =head1 DESCRIPTION
  482 
  483 This Perl module calulates percentage points (6 significant digits) of the u (standard normal) distribution,
  484 the student's t distribution, the chi-square distribution and the F distribution.
  485 
  486 It can also calculate the upper probability (6 significant digits) of the u (standard normal),
  487 the chi-square, the t and the F distribution.
  488 
  489 These critical values are needed to perform statistical tests, like the u test, the t test, the chi-squared test,
  490 and the F test, and to calculate confidence intervals.
  491 
  492 If you are interested in more precise algorithms you could look at:
  493  StatLib: http://lib.stat.cmu.edu/apstat/
  494  Applied Statistics Algorithms by Griffiths, P. and Hill, I.D., Ellis Horwood: Chichester (1985)
  495 
  496 =head1 AUTHOR
  497 
  498 Michael Kospach, mike.perl@gmx.at
  499 Nice formating, simplification and bug repair by Matthias Trautner Kromann, mtk@id.cbs.dk
  500 
  501 =cut
  502 
  503 
  504 
  505 
  506 
  507 
  508 
  509 
  510 

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9