program binom(binomp, xyin, xyplop, output);
(* produce the binom probabilities for a found black to white ratio
   Thomas Schneider *)

label 1; (* end of program *)

const
(* begin module version *)
version = 1.49; (* of binom.p 1995 Nov 12
origin 1995 nov 12 from binomial.p *)
(* end module version *)

(* begin module describe.binom *)
(*
name
   binom: produce the binomial probabilities for a found black to white ratio

synopsis
   binom(xyin: out, xyplop: out, binomp: in, output: out)

files

   binomp: parameters to control the program, on three lines:
      blacks and whites: two integers on the first line, representing the
         number of black balls and white balls obtained in an experiment
      probability of black, P
   output: messages to the user

description
   Suppose there exists a large bin containing both black and white
   balls.  The true probability of black balls in the bin is P, and
   the probatility of white balls is (1-P).  We obtain a sample of
   black and white balls from the bin, given as the first two parameters
   in binomp.  The probability of getting this black:white sample is:

                          (black+white)!  black      white
       p(black:white|P) = -------------- P      (1-P)
                          black!white!

   The program generates these probabilities for a given P.
   The results are in a form that the xyplo program can use to plot.

see also
   xyplo.p, binplo.p

author
   Thomas Dana Schneider

bugs
   none known

*)
(* end module describe.binom *)

var
   xyin, (* the input file to xyplop *)
   xyplop, (* the control file for xyplop *)
   binomp: (* control parameters for binom *)
      text;

(* begin module halt *) 
procedure halt; 
(* stop the program.  the procedure performs a goto to the end of the   
   program.  you must have a label: 
      label 1;
   declared, and also the end of the program must have this label:  
      1: end. 
   examples are in the module libraries.
   this is the only goto in the delila system. *) 
begin 
      writeln(output,' program halt.'); 
      goto 1
end;  
(* end module halt version = 'delmod 6.16 84 mar 12 tds/gds'; *)

(* begin module pbinomial *)
function pbinomial(black, white: integer; P: real): real;
(* Calculate the probability of obtaining a certain number of black and white
balls in a sample from a bin containing a given probability P of black balls.

This is the binomial distribution:
                       (black+white)!  black      white
    p(black:white|P) = -------------- P      (1-P)
                       black!white!
This is calulated by taking the log of the whole thing, figuring
out the sum of the parts and exponentiating. *)
var
   sum: real; (* log of the solution *)
   total: integer; (* sum of black and white *)

function combination(black, white: integer): real;
(* calculate the combinations for a number of black and white balls.
give the result as the log of the power *)
var index: integer; (* counting index *)
    sum: real; (* the sum of the logs of the factorial components *)
begin (* combination *)
   sum := 0;
   (* calculate the upper part of the factorial division.  This starts
      at one more than black and goes up to the total *)
   for index := (black+1) to (black+white) do begin
      sum := sum + ln(index)
   end;
   (* now remove the factorial of white *)
   for index := 1 to white do begin
      sum := sum - ln(index)
   end;
   combination := sum
end; (* combination *)

function power(black, white: integer; P: real): real;
(* calculate the power portion of the binomial distribution.
give the result as the log of the power *)
begin (* power *)
   if P >= 1.0
   then power := 0.0
   else power := black * ln(P) + white * ln(1-P)
end; (* power *)

begin
   total := black + white;
   sum := combination(black,white) + power(black,white,P);
   (* very small sum is essentially zero.  This avoids exponentiation
      of small values (which bombs). *)
   if sum < -50.0 then pbinomial := 0.0
                  else pbinomial := exp(sum);
end;
(* end module pbinomial *)

(* begin module tpbinom *)
procedure tpbinom (black, white: integer; P: real);
(* test the pbinomial procedure for the given P *)
var
   index: integer; (* counter for the running sum *)
   sum: real; (* the running sum of the probabilities *)
   total: integer; (* sum of black and white *)
begin
   writeln(output,'testing pbinomial procedure');
   sum := 0.0;
   total := black + white;
   for index := 0 to total do begin
      sum := sum + pbinomial(index,total-index,P);
   end;
   writeln(output,'sum of p(b,w|P) for ',total:1,' balls is ',sum:10:8)
end;
(* end module tpbinom *)

(* begin module binom.makexyplop *)
procedure makexyplop(var x: text; maxblack: integer);
(* write the xyplop values out *)
begin
writeln(x,'1 2       zerox zeroy');
writeln(x,'x 0 ',maxblack:1,'     zx min max          ');
writeln(x,'y 0 1 zy min max          ');
writeln(x,'20 10 1  1  xinterval yinterval ');
writeln(x,'4 5       xwidth    ywidth    ');
writeln(x,'0 3       xdecimal  ydecimal  ');
writeln(x,'6 6     xsize     ysize     ');
writeln(x,'number of blacks'); (* label x *)
writeln(x,'p(number of blacks)'); (* label y *)
(*writeln(x,'p(b=',black:1,': w=',white:1, ' | number b)'); old label *)
writeln(x,'a         no cross hairs');
writeln(x,'n         no log');
writeln(x,'n         no log');
writeln(x,' **********************************');
writeln(x,'1 2       xcolumn   ycolumn');
writeln(x,'0         symbol column');
writeln(x,'0 0       xscolumn   yscolumn');
writeln(x,'0 0 0     hue saturation brightness columns for color manipulation');
writeln(x,' **********************************');
writeln(x,'+         symbol-to-plot');
writeln(x,'+         symbol flag');
writeln(x,'0.01      symbol x size');
writeln(x,'0.01      symbol y size');
writeln(x,'cl 0.05   connected line  ');
writeln(x,'n 0.05    no regression-line  ');
writeln(x,' **********************************');
writeln(x,'.');
writeln(x,' **********************************');
end;
(* end module binom.makexyplop *)

(* begin module binom.themain *)
procedure themain(var xyin, xyplop, binomp: text);
(* the main procedure of the binom program *)
var
   black, white: integer; (* number of black and white balls found *)
   P: real; (* the probability of black balls *)
   Q: real; (* the probability of white balls *)
   mean: real; (* the mean of the distribution *)
   pbin: real; (* the calculated probability *)
   sigma: real; (* the standard deviation of the distribution *)
   total: integer; (* the total of blacks and whites *)
begin
   writeln(output,'binom ',version:4:2);

   reset(binomp);
   readln(binomp,black,white);
   readln(binomp,P);

   if (P < 0) then begin
      writeln(output, 'probability P must be positive');
      halt
   end;

   if (P > 1) then begin
      writeln(output, 'probability P must be less than 1.0');
      halt
   end;

   total := black + white;
   Q := 1-P;
   mean := total*P;
   sigma := sqrt(total*P*Q);

(*
   rewrite(xyplop);
*)

   rewrite(xyin);

   writeln(xyin,'* binom ',version:4:2);
   writeln(xyin,'* black=',black:1,': white=',white:1);
   writeln(xyin,'* probability of black =',P:5:3);
   writeln(xyin,'* mean  = ',mean:4:2);
   writeln(xyin,'* sigma = ',sigma:4:2);
   writeln(xyin,'* Zblack = ',((black - mean)/sigma):4:2);

(*
   writeln(xyin,'                  mean = ',mean:4:2);
   writeln(xyin,'                 sigma = ',sigma:4:2);
   writeln(xyin,'sharpness = mean/sigma = ',mean/sigma:4:1);
*)

   pbin := pbinomial(black,white,P);

   writeln(xyin, 'probability of this black/white combination: ',pbin:10:8)
end;
(* end module binom.themain *)

begin (* binom *)
   themain(xyin, xyplop, binomp);
1: end.