program binplo(xyin, xyplop, binplop, output);
(* produce the binomial probabilities for a found black to white ratio
   Thomas Dana Schneider  1987 *)

label 1; (* end of program *)

const
(* begin module version *)
version = 1.29; (* of binplo.p, 1987 feb 10 *)
(* end module version *)

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

synopsis
   binplo(xyin: out, xyplop: out, binplop: in, output: out)

files
   xyin: a table of probabilities of finding the given black to white
      ratio, versus the true probability.  The form is a series of lines
      that begin with '* ', followed by two columns of numbers.
      The first column is the value of fraction, and the second column is
      the corresponding value of p(black:white|fraction) = the probability of
      obtaining black and white given fraction.  This file is direct input
      to the xyplo program.
   xyplop: the controls for the xyplo program to generate the graph.
      These may be modified by the user before plotting.
   binplop: parameters to control the program
      blacks and whites: two integers on the first line, representing the
         number of black balls and white balls obtained in an experiment
      points: one integer on the second line, how many data points should
         be generated in the fout.  If points is zero, then the program
         tests its binomial probability procedure by adding all the
         probabilities that correspond to the binomial distribution.
         For example, with 1 black and 18 white balls, the test is to
         add the probabilities for (0,19), (1,18), ... (19,0).  This
         value should be close to 1.00 if the procedure is correct.

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

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

   the program generates these probabilities for all values of fraction,
   and gives the results in a form that the xyplo program can use to plot.

see also
   xyplo.p

author
   Thomas Dana Schneider

bugs
   none known

*)
(* end module describe.binplo *)

var
   xyin, (* the input file to xyplop *)
   xyplop, (* the control file for xyplop *)
   binplop: (* control parameters for binplo *)
      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; fraction: real): real;
(* calculate the probability of obtaining a certain number of black and
white balls in a sample from a bin containing a given fraction of balls.
This is the binomial distribution:
                              (black+white)!         black             white
    p(black:white|fraction) = -------------- fraction      (1-fraction)
                              black!white!
This is calulated by taking the log of the whole thing, figuring
out the sum of the parts and exponentiating. *)
var
   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; fraction: real): real;
(* calculate the power portion of the binomial distribution.
give the result as the log of the power *)
begin (* power *)
   power := black * ln(fraction) + white * ln(1-fraction)
end; (* power *)

begin
   total := black + white;
   pbinomial := exp(combination(black,white) + power(black,white,fraction));
end;
(* end module pbinomial *)

(* begin module tpbinom *)
procedure tpbinom (black, white: integer; fraction: real);
(* test the pbinomial procedure for the given fraction *)
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,fraction);
   end;
   writeln(output,'sum of p(b,w|fraction) for ',total:1,' balls is ',sum:10:8)
end;
(* end module tpbinom *)

(* begin module binplo.makexyplop *)
procedure makexyplop(var x: text; black, white: integer);
(* write the xyplop values out *)
begin
writeln(x,'x 0 1     zx min max          ');
writeln(x,'y 0 ',pbinomial(black,white,black/(black+white)):10:8,
                   ' zy min max          ');
writeln(x,'20 20     xinterval yinterval ');
writeln(x,'4 5       xwidth    ywidth    ');
writeln(x,'2 3       xdecimal  ydecimal  ');
writeln(x,'5.5 4     xsize     ysize     ');
writeln(x,'fraction of black'); (* label x *)
writeln(x,'p(b=',black:1,': w=',white:1,
          ' | fraction b)'); (* label y *)
writeln(x,'n         no cross hairs');
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,' --');
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 binplo.makexyplop *)

(* begin module binplo.themain *)
procedure themain(var xyin, xyplop, binplop: text);
(* the main procedure of the binplo program *)
var
   black, white: integer; (* number of black and white balls found *)
   fraction: real; (* the true fraction of black and white balls *)
   index: integer; (* index to the plotting of the points *)
   interval: real; (* 1/index *)
   points: integer; (* number of points to plot *)
   p: real; (* the calculated probability *)
begin
   writeln(output,'binplo ',version:4:2);

   reset(binplop);
   readln(binplop,black,white);
   readln(binplop,points);

   rewrite(xyplop);
   makexyplop(xyplop,black,white);

   rewrite(xyin);
   writeln(xyin,'* binplo ',version:4:2);
   writeln(xyin,'* ',points:1,' intervals plotted');
   writeln(xyin,'* black=',black:1,': white=',white:1);

   if points <= 0 then begin
      if points = 0 then tpbinom(black,white,black/(black+white))
      else begin
         writeln(output,' you cannot have a negative number of points');
         halt
      end
   end
   else begin
      interval := 1/points;
      fraction := 0.0;
      for index := 1 to points do begin
         fraction := fraction+interval;
         p := pbinomial(black,white,fraction);
         if trunc(p*1e8) > 0 then writeln(xyin,fraction:10:8,' ',p:10:8)
      end
   end
end;
(* end module binplo.themain *)

begin (* binplo *)
   themain(xyin, xyplop, binplop);
1: end.