#!/usr/bin/perl
# program greg5.prl 
# simplified version of greg3.prl and greg4.prl to print out listing
# that can be read into a stats package
# calculates probabilities of a maximum number of repeats when drawing
# n of t equiprobable objects with replacement
#
# copyright claimed by Chris Evans on 5.vii.96
# although I claim the copyright on this horrible bit of programming
# I am happy for anyone to use it however they like provided:
# a)   you don't make a profit out of it (inconceivable!)
# and
# b)   you make a reasonable attempt to tell me how to improve it,
#      do the same thing analytically, and what you've done with it
# and
# c)   you make reasonable acknowledgement of having used it if 
#      publishing something done using it.
#
# I can be reached by Email: C.Evans@sghms.ac.uk
#   'phone and 'fax: [+44|0] 181-725 2540
#   snailmail:  Psychotherapy Section, Dept. Gen. Psychiatry,
#               St. George's Hospital Medical School,
#               Cranmer Terrace, London SW17 0RE Britain
#

$t = 28;
$n = 11;

$curr = "1";
$freq_curr{$curr} = $t;  # there are always $t ways of selecting one object
$prob = 1/$t;
print "1 1 1\n";


for ($count = 2; $count <= $n; $count++) {
   undef(%score);
   $big = $t**$count;
   foreach $curr (keys (%freq_curr)) {
      &expand_curr($curr, $freq_curr{$curr});
      &sort_new(%new);
   }
   foreach $new (keys (%new)) {
      @curr = &split_curr($new);
      $score = $curr[0];
      $score{$score} += $new{$new};
   }
   %freq_curr = %new;
   undef(%new);
   &cum_score;
   &check_total;
}



exit;

# ==== formats ===============
format STDOUT_TOP = 
 score         n            proportion
------   ------------     ------------
.

format STDOUT =
@#####   @###########     @#.#########
$val, $score{$val}, $propn
.


format CUMOUT_TOP = 
... and for scoring at or above cutting score
   cut
 score         n            proportion
------   ------------     ------------
.

format CUMOUT =
@#####   @###########     @#.#########
$r, $cum_n, $propn
.




# ==== subroutines ===========



sub check_total {
   local ($total) = 0;
   foreach $val (keys %score) {
      $total += $score{$val};
   }
   if ($total != $big) {
      print "\nSomething wrong somewhere\n";
      print "total = $total <> expected total = $big\n";
      print "might be O.K., just that numeric limit for Perl on this machine exceeded";
   }
}


sub cum_score {
   for($r = $count; $r > 0; $r--) {
      $cum_n = 0;
      for ($s = $r; $s <= $n; $s++) {
         $cum_n += $score{$s};
      }
      $propn = $cum_n/$big;
      print "$count $r $propn\n";
   }
}


sub tot_score {
   local($new_score);
   # work out the score of the current selection
   $new_score = $curr[0];
   # increment the counter for that score
   $score{$new_score}++;
}   


sub split_curr {
   # takes a string, splits it on spaces, returns a list
   local($curr) = @_; 
   split(' ',$curr);
}

sub join_curr {
   # takes a list, joins it using spaces, returns a string
   local(@curr) = @_;
   local($string);
   $string = join(' ',@curr);
}


sub expand_curr {
   # takes a string and its frequency and works out what will be 
   # generated from that by selecting another object
   local($curr,$freq) = @_;
   local(@curr) = &split_curr($curr);
   $total = $freq*$t; # the number of scores that will be generated
   for ($index = 0; $index <= $#curr; $index++) {
      if ($total) {
         @new = @curr;
         $new[$index] ++;
         $new = &join_curr(@new);
         $new{$new} += $freq;
         $total -= $freq;
      }
   }
   if ($total) {
      @new = @curr;
      push(@new,'1');
      $new = &join_curr(@new);
      $new{$new} += $total;
   }
}


sub sort_new {
   local(%tmp);
   local ($curr);
   foreach $curr (keys (%new)) {
      local(@newcurr) = &split_curr($curr);
      @newcurr = reverse(sort(@newcurr));
      local($newcurr) = &join_curr(@newcurr);
      $tmp{$newcurr} += $new{$curr};
   }
   %new = %tmp;
}

# === text block ==========