modules/Math-Random/example.pl

   1 #!/usr/bin/perl -w
   2 use Math::Random qw(:all);
   3 use Carp;
   4 use strict;
   5
   6 my $max_choice = 17;
   7 my $ans = "";
   8 my $input = "";
   9 my @args = ();
  10 my @result = ();
  11
  12  TEST: while (1) {
  13      print " Enter number corresponding to choice:\n",
  14      "      (0) Exit this program\n",
  15      "      (1) Generate Chi-Square deviates\n",
  16      "      (2) Generate noncentral Chi-Square deviates\n",
  17      "      (3) Generate F deviates\n",
  18      "      (4) Generate noncentral F  deviates\n",
  19      "      (5) Generate random permutation\n",
  20      "      (6) Generate uniform integers\n",
  21      "      (7) Generate uniform reals\n",
  22      "      (8) Generate beta deviates\n",
  23      "      (9) Generate binomial outcomes\n",
  24      "     (10) Generate Poisson outcomes\n",
  25      "     (11) Generate exponential deviates\n",
  26      "     (12) Generate gamma deviates\n",
  27      "     (13) Generate multinomial outcomes\n",
  28      "     (14) Generate normal deviates\n",
  29      "     (15) Generate negative binomial outcomes\n",
  30      "     (16) Generate multivariate normal deviates\n",
  31      "     (17) Generate random permuted index\n";
  32      $ans = <stdin>;
  33      chomp $ans;
  34      $ans = int($ans);
  35      last TEST if $ans == 0;
  36      unless ($ans > 0 and $ans <= $max_choice) {
  37          print "Try one of (1, ... $max_choice)\n";
  38          next TEST;
  39      }
  40      print "Enter phrase to initialize seeds:\n";
  41      my $phrase = <stdin>;
  42      chomp $phrase;
  43      random_set_seed_from_phrase($phrase);
  44
  45      if ($ans == 1) {
  46          print "Enter (space-separated) N, DF for chi-square deviates:\n";
  47          $input = <stdin>;
  48          chomp $input;
  49          @args = split(/\s+/,$input);
  50          @result = random_chi_square(@args);
  51          shift(@args);
  52          print_results(sampstat(@result),scalar(@result),
  53                        true_stats('chis',@args));
  54          next TEST;
  55      }
  56
  57      if ($ans == 2) {
  58          print "Enter (space-separated) N, DF, NONC for ",
  59          "noncentral chi-square deviates:\n";
  60          $input = <stdin>;
  61          chomp $input;
  62          @args = split(/\s+/,$input);
  63          @result = random_noncentral_chi_square(@args);
  64          shift(@args);
  65          print_results(sampstat(@result),scalar(@result),
  66                        true_stats('ncch',@args));
  67          next TEST;
  68      }
  69
  70      if ($ans == 3) {
  71          print "Enter (space-separated) N, DFN, DFD for F deviates:\n";
  72          $input = <stdin>;
  73          chomp $input;
  74          @args = split(/\s+/,$input);
  75          @result = random_f(@args);
  76          shift(@args);
  77          print_results(sampstat(@result),scalar(@result),
  78                        true_stats('f',@args));
  79          next TEST;
  80      }
  81
  82      if ($ans == 4) {
  83          print "Enter (space-separated) N, DFN, DFD, NONC for ",
  84          "non-central F deviates:\n";
  85          $input = <stdin>;
  86          chomp $input;
  87          @args = split(/\s+/,$input);
  88          @result = random_noncentral_f(@args);
  89          shift(@args);
  90          print_results(sampstat(@result),scalar(@result),
  91                        true_stats('ncf',@args));
  92          next TEST;
  93      }
  94
  95      if ($ans == 5) {
  96          print "Enter (space-separated) list for random permutation:\n";
  97          $input = <stdin>;
  98          chomp $input;
  99          @args = split(/\s+/,$input);
 100          print "Result (' : ' separated):\n",
 101          join(" : ", random_permutation(@args)),"\n";
 102          next TEST;
 103      }
 104
 105      if ($ans == 6) {
 106          print "Enter (space-separated) Maximum integer, Replications ",
 107          "per integer:\n";
 108          $input = <stdin>;
 109          chomp $input;
 110          @args = split(/\s+/,$input);
 111          @result = random_uniform_integer($args[0] * $args[1], 1, $args[0]);
 112          my @totals = (0) x ($args[0] + 1);
 113          my $val = 0;
 114          foreach $val (@result) { $totals[$val]++; }
 115          shift @totals;
 116          print "Result (' : ' separated):\n",join(" : ", @totals),"\n";
 117          next TEST;
 118      }
 119
 120      if ($ans == 7) {
 121          print "Enter (space-separated) N, LOWER, UPPER for ",
 122          "uniform real deviates:\n";
 123          $input = <stdin>;
 124          chomp $input;
 125          @args = split(/\s+/,$input);
 126          @result = random_uniform(@args);
 127          shift(@args);
 128          $args[0] = 0 unless defined($args[0]);
 129          $args[1] = 1 unless defined($args[1]);
 130          print_results(sampstat(@result),scalar(@result),
 131                        true_stats('unif',@args));
 132          next TEST;
 133      }
 134
 135      if ($ans == 8) {
 136          print "Enter (space-separated) N, A, B for ",
 137          "beta deviates:\n";
 138          $input = <stdin>;
 139          chomp $input;
 140          @args = split(/\s+/,$input);
 141          @result = random_beta(@args);
 142          shift(@args);
 143          print_results(sampstat(@result),scalar(@result),
 144                        true_stats('beta',@args));
 145          next TEST;
 146      }
 147
 148      if ($ans == 9) {
 149          print "Enter (space-separated) N, NTrials, P for ",
 150          "binomial outcomes:\n";
 151          $input = <stdin>;
 152          chomp $input;
 153          @args = split(/\s+/,$input);
 154          @result = random_binomial(@args);
 155          shift(@args);
 156          print_results(sampstat(@result),scalar(@result),
 157                        true_stats('bin',@args));
 158          next TEST;
 159      }
 160
 161      if ($ans == 10) {
 162          print "Enter (space-separated) N, MU for ",
 163          "poisson outcomes:\n";
 164          $input = <stdin>;
 165          chomp $input;
 166          @args = split(/\s+/,$input);
 167          @result = random_poisson(@args);
 168          shift(@args);
 169          print_results(sampstat(@result),scalar(@result),
 170                        true_stats('pois',@args));
 171          next TEST;
 172      }
 173
 174      if ($ans == 11) {
 175          print "Enter (space-separated) N, AV for ",
 176          "exponential deviates:\n";
 177          $input = <stdin>;
 178          chomp $input;
 179          @args = split(/\s+/,$input);
 180          @result = random_exponential(@args);
 181          shift(@args);
 182          $args[0] = 1 unless defined($args[0]);
 183          print_results(sampstat(@result),scalar(@result),
 184                        true_stats('expo',@args));
 185          next TEST;
 186      }
 187
 188      if ($ans == 12) {
 189          print "Enter (space-separated) N, A, R for ",
 190          "gamma deviates:\n";
 191          $input = <stdin>;
 192          chomp $input;
 193          @args = split(/\s+/,$input);
 194          @result = random_gamma(@args);
 195          shift(@args);
 196          print_results(sampstat(@result),scalar(@result),
 197                        true_stats('gamm',@args));
 198          next TEST;
 199      }
 200
 201      if ($ans == 13) {
 202          print "Enter (space-separated) list of prob.s for categories ",
 203          "for multinomial outcomes:\n";
 204          my $input = <stdin>;
 205          chomp $input;
 206          @args = split(/\s+/,$input);
 207          print "Please enter number of events to be classified:\n";
 208          my $n = <stdin>;
 209          chomp $n;
 210          @result = random_multinomial($n,@args);
 211          print "Result:\n",join(" ", @result),"\n";
 212          my $sum = 0;
 213          my $val;
 214          foreach $val (@result) {$sum += $val; }
 215          foreach $val (@result) {$val /= $sum; }
 216          print "Observed proportions:\n",join(" ", @result),"\n";
 217          pop @args;
 218          $sum = 0;
 219          foreach $val (@args) {$sum += $val; }
 220          push @args, (1 - $sum);
 221          print "Expected proportions:\n",join(" ", @args),"\n";
 222          next TEST;
 223      }
 224
 225      if ($ans == 14) {
 226          print "Enter (space-separated) N, AV, SD for ",
 227          "normal deviates:\n";
 228          $input = <stdin>;
 229          chomp $input;
 230          @args = split(/\s+/,$input);
 231          @result = random_normal(@args);
 232          shift(@args);
 233          $args[0] = 0 unless defined($args[0]);
 234          $args[1] = 1 unless defined($args[1]);
 235          print_results(sampstat(@result),scalar(@result),
 236                        true_stats('norm',@args));
 237          next TEST;
 238      }
 239
 240      if ($ans == 15) {
 241          print "Enter (space-separated) N, NEvents, P for ",
 242          "negative binomial outcomes:\n";
 243          $input = <stdin>;
 244          chomp $input;
 245          @args = split(/\s+/,$input);
 246          @result = random_negative_binomial(@args);
 247          shift(@args);
 248          print_results(sampstat(@result),scalar(@result),
 249                        true_stats('nbin',@args));
 250          next TEST;
 251      }
 252
 253      if ($ans == 16) {
 254          print "Enter dimension of multivariate deviate:\n";
 255          my $p = <stdin>;
 256          chomp $p;
 257          print "Enter mean vector of length $p (space separated):\n";
 258          my $temp = <stdin>;
 259          chomp $temp;
 260          $temp =~ s/^\s*//;
 261          my @mean;
 262          @mean = split(/\s+/,$temp);
 263          print "Enter (symmetric, $p by $p) covariance matrix\n",
 264          "One space-separated row per line:\n";
 265          my $val;
 266          my @covariance = (0) x $p;
 267          foreach $val (@covariance) {
 268              $temp = <stdin>;
 269              chomp $temp;
 270              $temp =~ s/^\s*//;
 271              $val = [ split(/\s+/,$temp) ];
 272          }
 273          print "Enter number of observations:\n";
 274          my $n = <stdin>;
 275          chomp $n;
 276          my @ans = random_multivariate_normal($n,@mean,@covariance);
 277          my $template = join(" ",('%15.7g') x $p);
 278          print "\nResults:\n";
 279          foreach $val (@ans) {
 280              printf "$template\n",@{$val};
 281          }
 282      }
 283
 284      if ($ans == 17) {
 285          print "Enter N (size of array) for random permuted index:\n";
 286          $input = <stdin>;
 287          chomp $input;
 288          print "Result (' : ' separated):\n",
 289          join(" : ", random_permuted_index($input)),"\n";
 290          next TEST;
 291      }
 292  }
 293 print "Normal termination of tester.\n";
 294
 295 sub sampstat { # gets sample statistics for array - returns array of stats
 296     my $n = scalar(@_);
 297     return () unless $n > 0;
 298     return ($_[0], 0, $_[0], $_[0]) unless $n > 1;
 299     my($min) = my($max) = $_[0];
 300     my $val;
 301     my $avg = 0;
 302     foreach $val (@_) {
 303         $avg += $val;
 304         $min = $val if $val < $min;
 305         $max = $val if $val > $max;
 306     }
 307     $avg /= $n;
 308     my $var = 0;
 309     foreach $val (@_) {
 310         $var += ($val - $avg)**2;
 311     }
 312     $var /= ($n - 1);
 313     return ($avg, $var, $min, $max);
 314 }
 315
 316 sub print_results { # Arguments: $avg, $var, $min, $max, $nobs, $travg, $trvar
 317     my($avg, $var, $min, $max, $nobs, $travg, $trvar) = @_;
 318     print "Results:\n";
 319     printf "Number of observations: %d\n", $nobs;
 320     printf "Mean    : %15.7g   True    : %15.7g\n", $avg, $travg;
 321     printf "Variance: %15.7g   True    : %15.7g\n", $var, $trvar;
 322     printf "Minimum : %15.7g   Maximum : %15.7g\n", $min, $max;
 323 }
 324
 325 sub true_stats { # Arguments: $type, @parin; Returns: $av, $var
 326
 327 ########################################################################
 328 #     Returns mean and variance for a number of statistical distribution
 329 #     as a function of their parameters.
 330 #
 331 #
 332 #                              Arguments
 333 #
 334 #
 335 #    $type --> Character string indicating type of distribution
 336 #             'chis' chisquare
 337 #             'ncch' noncentral chisquare
 338 #             'f'    F (variance ratio)
 339 #             'ncf'  noncentral f
 340 #             'unif' uniform
 341 #             'beta' beta distribution
 342 #             'bin'  binomial
 343 #             'pois' poisson
 344 #             'expo' exponential
 345 #             'gamm' gamma
 346 #             'norm' normal
 347 #             'nbin' negative binomial
 348 #
 349 #    @parin --> Array containing parameters of distribution
 350 #              chisquare
 351 #               $parin[0] is df
 352 #              noncentral chisquare
 353 #               $parin[0] is df
 354 #               $parin[1] is noncentrality parameter
 355 #              F (variance ratio)
 356 #               $parin[0] is df numerator
 357 #               $parin[1] is df denominator
 358 #              noncentral F
 359 #               $parin[0] is df numerator
 360 #               $parin[1] is df denominator
 361 #               $parin[2] is noncentrality parameter
 362 #              uniform
 363 #               $parin[0] is LOW bound
 364 #               $parin[1] is HIGH bound
 365 #              beta
 366 #               $parin[0] is A
 367 #               $parin[1] is B
 368 #              binomial
 369 #               $parin[0] is Number of trials
 370 #               $parin[1] is Prob Event at Each Trial
 371 #              poisson
 372 #               $parin[0] is Mean
 373 #              exponential
 374 #               $parin[0] is Mean
 375 #              gamma
 376 #               $parin[0] is A
 377 #               $parin[1] is R
 378 #              normal
 379 #               $parin[0] is Mean
 380 #               $parin[1] is Standard Deviation
 381 #              negative binomial
 382 #               $parin[0] is required Number of events
 383 #               $parin[1] is Probability of event
 384 #
 385 #     $av <-- Mean of specified distribution with specified parameters
 386 #
 387 #     $var <-- Variance of specified distribution with specified paramete
 388 #
 389 #
 390 #                              Note
 391 #
 392 #
 393 #     $av and $var will be returned -1 if mean or variance is infinite
 394 #
 395 #**********************************************************************
 396     my $type = shift(@_);
 397     my @parin = @_;
 398     my($av, $var, $a, $b, $range) = (-1,-1,0,0,0);
 399
 400   TYPE: {
 401       if (('chis') eq ($type)){
 402           $av = $parin[0];
 403           $var = 2.0*$parin[0];
 404           last TYPE;}
 405
 406       if (('ncch') eq ($type)) {
 407           $a = $parin[0] + $parin[1];
 408           $b = $parin[1]/$a;
 409           $av = $a;
 410           $var = 2.0*$a* (1.0+$b);
 411           last TYPE;}
 412
 413       if (('f') eq ($type)) {
 414           unless ($parin[1] <= 2.0001) {
 415               $av = $parin[1]/ ($parin[1]-2.0);
 416           }
 417           unless ($parin[1] <= 4.0001) {
 418               $var = (2.0*$parin[1]**2* ($parin[0]+$parin[1]-2.0))/
 419                   ($parin[0]* ($parin[1]-2.0)**2* ($parin[1]-4.0));
 420           }
 421           last TYPE;}
 422
 423       if (('ncf') eq ($type)) {
 424           unless ($parin[1] <= 2.0001){
 425               $av = ($parin[1]* ($parin[0]+$parin[2]))/
 426                   (($parin[1]-2.0)*$parin[0]);
 427           }
 428           unless ($parin[1] <= 4.0001) {
 429               $a = ($parin[0]+$parin[2])**2 + ($parin[0]+2.0*$parin[2])*
 430                   ($parin[1]-2.0);
 431               $b = ($parin[1]-2.0)**2* ($parin[1]-4.0);
 432               $var = 2.0* ($parin[1]/$parin[0])**2* ($a/$b);
 433           }
 434           last TYPE;}
 435
 436       if (('unif') eq ($type)) {
 437           $range = $parin[1] - $parin[0];
 438           $av = $parin[0] + $range/2.0;
 439           $var = $range**2/12.0;
 440           last TYPE;}
 441
 442       if (('beta') eq ($type)) {
 443           $av = $parin[0]/ ($parin[0]+$parin[1]);
 444           $var = ($av*$parin[1])/ (($parin[0]+$parin[1])*
 445                                    ($parin[0]+$parin[1]+1.0));
 446           last TYPE;}
 447
 448       if (('bin') eq ($type)) {
 449           $av = $parin[0]*$parin[1];
 450           $var = $av* (1.0-$parin[1]);
 451           last TYPE;}
 452
 453       if (('pois') eq ($type)) {
 454           $av = $parin[0];
 455           $var = $parin[0];
 456           last TYPE;}
 457
 458       if (('expo') eq ($type)) {
 459           $av = $parin[0];
 460           $var = $parin[0]**2;
 461           last TYPE;}
 462
 463       if (('gamm') eq ($type)) {
 464           $av = $parin[1] / $parin[0];
 465           $var = $av / $parin[0];
 466           last TYPE;}
 467
 468       if (('norm') eq ($type)) {
 469           $av = $parin[0];
 470           $var = $parin[1]**2;
 471           last TYPE;}
 472
 473       if (('nbin') eq ($type)) {
 474           $av = $parin[0] * (1.0 - $parin[1]) / $parin[1];
 475           $var = $av / $parin[1];
 476           last TYPE;}
 477
 478       croak "Unimplemented \$type: $type in true_stats";
 479   }
 480     return ($av,$var);
 481 }