## Compute Adjusted Wald 95% CI for binary results
## Assumption is that success rate will usually be > 0.5


###  ###  ###  ###  ###  ###  ###  ###  ###  ###  ###
   sub getnum
   {
      my ($prompt) = @_;
      my ($inp);

      print $prompt;
      $inp = <STDIN>;
      if ($inp =~ m%(\d+)%)
      {  return $1; }
      else
      {  die "Bad input\n"; }
   }


###  ###  ###  ###  ###  ###  ###  ###  ###  ###  ###
sub adj_wald
{  # compute adjusted Wald CI

   my ($succ, $trials) = @_;
   my ($signif_95, $signif_90, $sig2, $base_rate, $p_center, $p_hat,
       $half_CI, $lo_CI, $hi_CI );

   $signif_95 = 1.96;   # set for 95% confidence
   $signif_90 = 1.645;  # set for 90% confidence
   $sig2   = $signif_95 * $signif_95;

   # estimate of successful ratio
   $base_rate = $succ / $trials;

   # Find CI center
   $p_center = ($succ + $sig2 / 2) / ($trials + $sig2);

   # Adjusted estimate
   if ($base_rate <= 0.5)
   {  $p_hat = $p_center; }
   elsif ($base_rate <= 0.9)
   {  $p_hat = $base_rate;  }
   else
   {  $p_hat = ($succ + 1) / ($trials + 2); }

##   printf RFILE "base = %8.4f; cent = %8.4f; hat = %8.4f; \n\n",
##      $base_rate, $p_center, $p_hat;

   if ($succ < $trials)
   {
      $half_CI = $signif_95 * sqrt ($p_center * (1 - $p_center) / ($trials + $sig2));

      $lo_CI = $p_center - $half_CI;
      $hi_CI = $p_center + $half_CI;
      $hi_CI = 1 if ($hi_CI > 1);
   }
   else
   {  #special handling for succ = trials

      $half_CI = $signif_90 * sqrt ($p_center * (1 - $p_center) / ($trials + $sig2));

      $lo_CI = $p_hat - $half_CI;
      $hi_CI = 1;
   }

   ## print RFILE "half_CI = $half_CI \n";
   return ($p_hat, $lo_CI, $hi_CI);

}

###  ###  ###  ###  ###  ###  ###  ###  ###  ###  ###
###  ###  ###   MAIN LINE STARTS HERE   ###  ###  ###
###  ###  ###  ###  ###  ###  ###  ###  ###  ###  ###

                                                                  
   print "\nCompute Adjusted Wald 95% CI for binary results \n\n";

   $succ   = getnum ("Enter # successes: ");
   $trials = getnum ("Enter Total # trials: ");

   if ($succ <= 0)
   {  die "#successes must be > 0 \n";
   }
   if ($trials <= 0)
   {  die "#trials must be > 0 \n";
   }
   if ($succ > $trials)
   {  die "Cannot handle #successes > #trials. \n";
   }

   ($p_hat, $lo_CI, $hi_CI) = adj_wald ($succ, $trials);

   printf "Estimated success rate = %8.4f\n\n95% CI = [ %8.4f, %8.4f ] \n",
      $p_hat, $lo_CI, $hi_CI;