#!/usr/bin/perl -w
use Math::Random qw(:all);
use Carp;
use strict;
my $max_choice = 17;
my $ans = "";
my $input = "";
my @args = ();
my @result = ();
TEST: while (1) {
print " Enter number corresponding to choice:\n",
" (0) Exit this program\n",
" (1) Generate Chi-Square deviates\n",
" (2) Generate noncentral Chi-Square deviates\n",
" (3) Generate F deviates\n",
" (4) Generate noncentral F deviates\n",
" (5) Generate random permutation\n",
" (6) Generate uniform integers\n",
" (7) Generate uniform reals\n",
" (8) Generate beta deviates\n",
" (9) Generate binomial outcomes\n",
" (10) Generate Poisson outcomes\n",
" (11) Generate exponential deviates\n",
" (12) Generate gamma deviates\n",
" (13) Generate multinomial outcomes\n",
" (14) Generate normal deviates\n",
" (15) Generate negative binomial outcomes\n",
" (16) Generate multivariate normal deviates\n",
" (17) Generate random permuted index\n";
$ans = <stdin>;
chomp $ans;
$ans = int($ans);
last TEST if $ans == 0;
unless ($ans > 0 and $ans <= $max_choice) {
print "Try one of (1, ... $max_choice)\n";
next TEST;
}
print "Enter phrase to initialize seeds:\n";
my $phrase = <stdin>;
chomp $phrase;
random_set_seed_from_phrase($phrase);
if ($ans == 1) {
print "Enter (space-separated) N, DF for chi-square deviates:\n";
$input = <stdin>;
chomp $input;
@args = split(/\s+/,$input);
@result = random_chi_square(@args);
shift(@args);
print_results(sampstat(@result),scalar(@result),
true_stats('chis',@args));
next TEST;
}
if ($ans == 2) {
print "Enter (space-separated) N, DF, NONC for ",
"noncentral chi-square deviates:\n";
$input = <stdin>;
chomp $input;
@args = split(/\s+/,$input);
@result = random_noncentral_chi_square(@args);
shift(@args);
print_results(sampstat(@result),scalar(@result),
true_stats('ncch',@args));
next TEST;
}
if ($ans == 3) {
print "Enter (space-separated) N, DFN, DFD for F deviates:\n";
$input = <stdin>;
chomp $input;
@args = split(/\s+/,$input);
@result = random_f(@args);
shift(@args);
print_results(sampstat(@result),scalar(@result),
true_stats('f',@args));
next TEST;
}
if ($ans == 4) {
print "Enter (space-separated) N, DFN, DFD, NONC for ",
"non-central F deviates:\n";
$input = <stdin>;
chomp $input;
@args = split(/\s+/,$input);
@result = random_noncentral_f(@args);
shift(@args);
print_results(sampstat(@result),scalar(@result),
true_stats('ncf',@args));
next TEST;
}
if ($ans == 5) {
print "Enter (space-separated) list for random permutation:\n";
$input = <stdin>;
chomp $input;
@args = split(/\s+/,$input);
print "Result (' : ' separated):\n",
join(" : ", random_permutation(@args)),"\n";
next TEST;
}
if ($ans == 6) {
print "Enter (space-separated) Maximum integer, Replications ",
"per integer:\n";
$input = <stdin>;
chomp $input;
@args = split(/\s+/,$input);
@result = random_uniform_integer($args[0] * $args[1], 1, $args[0]);
my @totals = (0) x ($args[0] + 1);
my $val = 0;
foreach $val (@result) { $totals[$val]++; }
shift @totals;
print "Result (' : ' separated):\n",join(" : ", @totals),"\n";
next TEST;
}
if ($ans == 7) {
print "Enter (space-separated) N, LOWER, UPPER for ",
"uniform real deviates:\n";
$input = <stdin>;
chomp $input;
@args = split(/\s+/,$input);
@result = random_uniform(@args);
shift(@args);
$args[0] = 0 unless defined($args[0]);
$args[1] = 1 unless defined($args[1]);
print_results(sampstat(@result),scalar(@result),
true_stats('unif',@args));
next TEST;
}
if ($ans == 8) {
print "Enter (space-separated) N, A, B for ",
"beta deviates:\n";
$input = <stdin>;
chomp $input;
@args = split(/\s+/,$input);
@result = random_beta(@args);
shift(@args);
print_results(sampstat(@result),scalar(@result),
true_stats('beta',@args));
next TEST;
}
if ($ans == 9) {
print "Enter (space-separated) N, NTrials, P for ",
"binomial outcomes:\n";
$input = <stdin>;
chomp $input;
@args = split(/\s+/,$input);
@result = random_binomial(@args);
shift(@args);
print_results(sampstat(@result),scalar(@result),
true_stats('bin',@args));
next TEST;
}
if ($ans == 10) {
print "Enter (space-separated) N, MU for ",
"poisson outcomes:\n";
$input = <stdin>;
chomp $input;
@args = split(/\s+/,$input);
@result = random_poisson(@args);
shift(@args);
print_results(sampstat(@result),scalar(@result),
true_stats('pois',@args));
next TEST;
}
if ($ans == 11) {
print "Enter (space-separated) N, AV for ",
"exponential deviates:\n";
$input = <stdin>;
chomp $input;
@args = split(/\s+/,$input);
@result = random_exponential(@args);
shift(@args);
$args[0] = 1 unless defined($args[0]);
print_results(sampstat(@result),scalar(@result),
true_stats('expo',@args));
next TEST;
}
if ($ans == 12) {
print "Enter (space-separated) N, A, R for ",
"gamma deviates:\n";
$input = <stdin>;
chomp $input;
@args = split(/\s+/,$input);
@result = random_gamma(@args);
shift(@args);
print_results(sampstat(@result),scalar(@result),
true_stats('gamm',@args));
next TEST;
}
if ($ans == 13) {
print "Enter (space-separated) list of prob.s for categories ",
"for multinomial outcomes:\n";
my $input = <stdin>;
chomp $input;
@args = split(/\s+/,$input);
print "Please enter number of events to be classified:\n";
my $n = <stdin>;
chomp $n;
@result = random_multinomial($n,@args);
print "Result:\n",join(" ", @result),"\n";
my $sum = 0;
my $val;
foreach $val (@result) {$sum += $val; }
foreach $val (@result) {$val /= $sum; }
print "Observed proportions:\n",join(" ", @result),"\n";
pop @args;
$sum = 0;
foreach $val (@args) {$sum += $val; }
push @args, (1 - $sum);
print "Expected proportions:\n",join(" ", @args),"\n";
next TEST;
}
if ($ans == 14) {
print "Enter (space-separated) N, AV, SD for ",
"normal deviates:\n";
$input = <stdin>;
chomp $input;
@args = split(/\s+/,$input);
@result = random_normal(@args);
shift(@args);
$args[0] = 0 unless defined($args[0]);
$args[1] = 1 unless defined($args[1]);
print_results(sampstat(@result),scalar(@result),
true_stats('norm',@args));
next TEST;
}
if ($ans == 15) {
print "Enter (space-separated) N, NEvents, P for ",
"negative binomial outcomes:\n";
$input = <stdin>;
chomp $input;
@args = split(/\s+/,$input);
@result = random_negative_binomial(@args);
shift(@args);
print_results(sampstat(@result),scalar(@result),
true_stats('nbin',@args));
next TEST;
}
if ($ans == 16) {
print "Enter dimension of multivariate deviate:\n";
my $p = <stdin>;
chomp $p;
print "Enter mean vector of length $p (space separated):\n";
my $temp = <stdin>;
chomp $temp;
$temp =~ s/^\s*//;
my @mean;
@mean = split(/\s+/,$temp);
print "Enter (symmetric, $p by $p) covariance matrix\n",
"One space-separated row per line:\n";
my $val;
my @covariance = (0) x $p;
foreach $val (@covariance) {
$temp = <stdin>;
chomp $temp;
$temp =~ s/^\s*//;
$val = [ split(/\s+/,$temp) ];
}
print "Enter number of observations:\n";
my $n = <stdin>;
chomp $n;
my @ans = random_multivariate_normal($n,@mean,@covariance);
my $template = join(" ",('%15.7g') x $p);
print "\nResults:\n";
foreach $val (@ans) {
printf "$template\n",@{$val};
}
}
if ($ans == 17) {
print "Enter N (size of array) for random permuted index:\n";
$input = <stdin>;
chomp $input;
print "Result (' : ' separated):\n",
join(" : ", random_permuted_index($input)),"\n";
next TEST;
}
}
print "Normal termination of tester.\n";
sub sampstat { # gets sample statistics for array - returns array of stats
my $n = scalar(@_);
return () unless $n > 0;
return ($_[0], 0, $_[0], $_[0]) unless $n > 1;
my($min) = my($max) = $_[0];
my $val;
my $avg = 0;
foreach $val (@_) {
$avg += $val;
$min = $val if $val < $min;
$max = $val if $val > $max;
}
$avg /= $n;
my $var = 0;
foreach $val (@_) {
$var += ($val - $avg)**2;
}
$var /= ($n - 1);
return ($avg, $var, $min, $max);
}
sub print_results { # Arguments: $avg, $var, $min, $max, $nobs, $travg, $trvar
my($avg, $var, $min, $max, $nobs, $travg, $trvar) = @_;
print "Results:\n";
printf "Number of observations: %d\n", $nobs;
printf "Mean : %15.7g True : %15.7g\n", $avg, $travg;
printf "Variance: %15.7g True : %15.7g\n", $var, $trvar;
printf "Minimum : %15.7g Maximum : %15.7g\n", $min, $max;
}
sub true_stats { # Arguments: $type, @parin; Returns: $av, $var
########################################################################
# Returns mean and variance for a number of statistical distribution
# as a function of their parameters.
#
#
# Arguments
#
#
# $type --> Character string indicating type of distribution
# 'chis' chisquare
# 'ncch' noncentral chisquare
# 'f' F (variance ratio)
# 'ncf' noncentral f
# 'unif' uniform
# 'beta' beta distribution
# 'bin' binomial
# 'pois' poisson
# 'expo' exponential
# 'gamm' gamma
# 'norm' normal
# 'nbin' negative binomial
#
# @parin --> Array containing parameters of distribution
# chisquare
# $parin[0] is df
# noncentral chisquare
# $parin[0] is df
# $parin[1] is noncentrality parameter
# F (variance ratio)
# $parin[0] is df numerator
# $parin[1] is df denominator
# noncentral F
# $parin[0] is df numerator
# $parin[1] is df denominator
# $parin[2] is noncentrality parameter
# uniform
# $parin[0] is LOW bound
# $parin[1] is HIGH bound
# beta
# $parin[0] is A
# $parin[1] is B
# binomial
# $parin[0] is Number of trials
# $parin[1] is Prob Event at Each Trial
# poisson
# $parin[0] is Mean
# exponential
# $parin[0] is Mean
# gamma
# $parin[0] is A
# $parin[1] is R
# normal
# $parin[0] is Mean
# $parin[1] is Standard Deviation
# negative binomial
# $parin[0] is required Number of events
# $parin[1] is Probability of event
#
# $av <-- Mean of specified distribution with specified parameters
#
# $var <-- Variance of specified distribution with specified paramete
#
#
# Note
#
#
# $av and $var will be returned -1 if mean or variance is infinite
#
#**********************************************************************
my $type = shift(@_);
my @parin = @_;
my($av, $var, $a, $b, $range) = (-1,-1,0,0,0);
TYPE: {
if (('chis') eq ($type)){
$av = $parin[0];
$var = 2.0*$parin[0];
last TYPE;}
if (('ncch') eq ($type)) {
$a = $parin[0] + $parin[1];
$b = $parin[1]/$a;
$av = $a;
$var = 2.0*$a* (1.0+$b);
last TYPE;}
if (('f') eq ($type)) {
unless ($parin[1] <= 2.0001) {
$av = $parin[1]/ ($parin[1]-2.0);
}
unless ($parin[1] <= 4.0001) {
$var = (2.0*$parin[1]**2* ($parin[0]+$parin[1]-2.0))/
($parin[0]* ($parin[1]-2.0)**2* ($parin[1]-4.0));
}
last TYPE;}
if (('ncf') eq ($type)) {
unless ($parin[1] <= 2.0001){
$av = ($parin[1]* ($parin[0]+$parin[2]))/
(($parin[1]-2.0)*$parin[0]);
}
unless ($parin[1] <= 4.0001) {
$a = ($parin[0]+$parin[2])**2 + ($parin[0]+2.0*$parin[2])*
($parin[1]-2.0);
$b = ($parin[1]-2.0)**2* ($parin[1]-4.0);
$var = 2.0* ($parin[1]/$parin[0])**2* ($a/$b);
}
last TYPE;}
if (('unif') eq ($type)) {
$range = $parin[1] - $parin[0];
$av = $parin[0] + $range/2.0;
$var = $range**2/12.0;
last TYPE;}
if (('beta') eq ($type)) {
$av = $parin[0]/ ($parin[0]+$parin[1]);
$var = ($av*$parin[1])/ (($parin[0]+$parin[1])*
($parin[0]+$parin[1]+1.0));
last TYPE;}
if (('bin') eq ($type)) {
$av = $parin[0]*$parin[1];
$var = $av* (1.0-$parin[1]);
last TYPE;}
if (('pois') eq ($type)) {
$av = $parin[0];
$var = $parin[0];
last TYPE;}
if (('expo') eq ($type)) {
$av = $parin[0];
$var = $parin[0]**2;
last TYPE;}
if (('gamm') eq ($type)) {
$av = $parin[1] / $parin[0];
$var = $av / $parin[0];
last TYPE;}
if (('norm') eq ($type)) {
$av = $parin[0];
$var = $parin[1]**2;
last TYPE;}
if (('nbin') eq ($type)) {
$av = $parin[0] * (1.0 - $parin[1]) / $parin[1];
$var = $av / $parin[1];
last TYPE;}
croak "Unimplemented \$type: $type in true_stats";
}
return ($av,$var);
}