#!/usr/bin/perl
# program greg.prl 
# calculates probabilities of a maximum number of repeats when drawing
# n of t equiprobable objects with replacement
#
# copyright claimed by Chris Evans on 3.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
#

# set the parameters
$t = 28;
$n = 4;

print $head_text;     # this explains what's going on, see end of program

print "\n\n  n = $n,  t = $t";
$big = $t**$n;        # total number of possible selections
	
undef @curr;
for ($r = 1; $r <= $n; $r++) {
   push(@curr,1);
}

for ($r = 1; $r <= $n; $r++) {
   $score{$r} = 0;
}


$which = 0;
&increment($which);

print "\n\nFinal score:\n";
&print_score;
&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 "\nO.K. total = $total = expected total = $big\n\n";
   } else {
      print "\nSomething wrong somewhere\n";
      print "\nO.K. total = $total <> expected total = $big\n\n";
   }
}


sub print_score {
   foreach $val (keys %score) {
      $propn = $score{$val}/$big;
      $perc  = $propn * 100;
      write;
   }
}

sub cum_score {
   print "\n\n\n";
   select((select(STDOUT), $| = 1, $- = 0, $^ = 'CUMOUT_TOP', $~ = 'CUMOUT')[0]);
   for($r = $n; $r > 0; $r--) {
      $cum_n = 0;
      for ($s = $r; $s <= $n; $s++) {
         $cum_n += $score{$s};
      }
      $propn = $cum_n/$big;
      write;
   }
   select((select(STDOUT), $| = 1, $- = 0, $^ = 'STDOUT_TOP', $~ = 'STDOUT')[0]);
}


sub tot_score {
   local($new_score);
   $new_score = &calc_score(@curr);
   $score{$new_score}++;
}   

sub calc_score {
  local(@curr) = @_;
  local(%n);
  local($r,$s,$max);
  for ($r = 0; $r < $n; $r++) {
     $s = $curr[$r];
     $n{$s}++;
  }
  undef $max;
  for ($r = 1; $r <= $t; $r++) {
     if ($n{$r} > $max) { $max = $n{$r} };
  }
  $max;
}


sub increment {
   local($which) = @_;
   local($next) = $which + 1;
   local($r);

   if ($which == ($n - 1)) { 
      for ($r = 1; $r <= $t; $r++) {
          $curr[$which] = $r;
          &tot_score;
      }
   } else {
      for ($r = 1; $r <= $t; $r++) {
         $curr[$which] = $r;
         &increment($next);
      }
   }
}




$head_text   = <<END;
Figures for highest number of identical objects in the
selection when n objects are selected from t equiprobable 
objects selecting with replacement.  E.g. if n = 8 and
there are 80 equipossible choices on each of the five
selections, say one object was selected three times,
another twice, another twice and three other objects
were selected once each then the score is 3.  Calculations
done by brute force simulation in Perl.  
(V. bad) programming by Chris Evans (C.Evans@sghms.ac.uk).
I would love someone to improve the programming and/or 
put me onto analytical solutions for this problem.
END