program binomial(xyin, xyplop, binomialp, output);
(* produce the binomial probabilities for a found black to white ratio

  Tom Schneider
  NCI/FCRDC Bldg 469. Room 144
  P.O. Box B
  Frederick, MD  21702-1201
  (301) 846-5581 (-5532 for messages)
  permanent email: toms@alum.mit.edu
  toms@ncifcrf.gov
  http://www-lecb.ncifcrf.gov/~toms/

  National Cancer Institute
  Laboratory of Experimental and Computational Biology

*)

label 1; (* end of program *)

const
(* begin module version *)
   version = 1.50; (* of binomial.p 2003 Sep 4
2003 Sep  4, 1.50:  upgrade to new xyplop parameters
1997 Dec 23, 1.49: previous version
1988 feb 24, 1.00: approximate origin
*)
(* end module version *)

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

synopsis
   binomial(xyin: out, xyplop: out, binomialp: in, output: out)

files

   xyin: a table of probabilities of finding the given white to black ratio,
      versus the probability.  The form is a series of comment lines that
      begin with '* ', followed by two columns of numbers.

      The first column is the number of whites.

      The second column is the corresponding value of p(white:black|pw) = the
      probability of obtaining the number of white and black given pw, the
      probability of white.  This file is direct input to the xyplo program.

      The third column is the cumulative sum of the second column.

   xyplop: the controls for the xyplo program to generate the graph.
      These may be modified by the user before plotting.

   binomialp: parameters to control the program, on three lines:

      balls: one integer on the first line, representing the total number of
	 black and white balls obtained in an experiment

      probability of white in the population.

      plot max: maximum number of balls to compute.  This allows one
         to limit the distribution plotted.

description

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

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

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

examples
   From Breiman1969 pages 24 and 25:

   The graphs on these pages can be simulated.

   Problem 6 a: "Compute the probability of getting exactly 10 Heads
   in 20 tosses of a fair coin. (Ans. .18)"

   For the binomialp file:

20      total number of black balls and white balls in an experiment
0.5     probability of drawing a white
100     plot max: maximum number of whites to show.

   Along with the rest of the distribution, program gives in xyin:
   10 0.1761970520020

documentation

@book{Breiman1969,
author = "L. Breiman",
title = "Probability and Stochastic Processes:
With a View Toward Applications",
publisher = "Houghton Mifflin Company",
address = "Boston",
year = "1969"}

see also
   xyplo.p, binplo.p

author
   Thomas Dana Schneider

bugs
   none known

*)
(* end module describe.binomial *)

var
   xyin, (* the input file to xyplop *)
   xyplop, (* the control file for xyplop *)
   binomialp: (* control parameters for binomial *)
      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
   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; fraction: real): real;
(* calculate the power portion of the binomial distribution.
give the result as the log of the power *)
begin (* power *)
   if fraction >= 1.0
   then power := 0.0
   else power := black * ln(fraction) + white * ln(1-fraction)
end; (* power *)

begin
   total := black + white;
   sum := combination(black,white) + power(black,white,fraction);
   (* 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; 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 binomial.makexyplop *)
procedure makexyplop(var x: text; maxwhite: integer);
(* write the xyplop values out *)
var
   xintervals: integer; (* for computing intervals to show *)
begin
   if maxwhite > 20
   then xintervals := 10
   else xintervals := maxwhite;
writeln(x,'3 6       zerox zeroy');
writeln(x,'zx 0 ',maxwhite:1,'     zx min max          ');
writeln(x,'zy 0 1 zy min max          ');
writeln(x,xintervals:1,' 10 1 1 xinterval yinterval ');
writeln(x,'6 6       xwidth    ywidth    ');
writeln(x,'0 3       xdecimal  ydecimal  ');
writeln(x,'15 15   xsize     ysize     ');
writeln(x,'number of white'); (* label x *)
writeln(x,'p(number of white)'); (* label y *)
(*writeln(x,'p(b=',black:1,': w=',white:1, ' | number b)'); old label *)
writeln(x,'a         no cross hairs');
writeln(x,'n 2       no log');
writeln(x,'n 2       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,'c         symbol-to-plot');
writeln(x,'o         symbol flag');
writeln(x,'0.250     symbol x size');
writeln(x,'0.250     symbol y size');
writeln(x,'cl 0.05   connected line  ');
writeln(x,'n 0.05    no regression-line  ');
writeln(x,' *-');
writeln(x,'.');
writeln(x,' *-');
writeln(x,' *-');
writeln(x,'p 1.50 0.50 1.50 1.00 edgecontrol (p=page),',
          ' edgeleft, edgeright, edgelow, edgehigh in cm');
writeln(x,'9.01 version of xyplo that this parameter file is designed for.');
end;
(* end module binomial.makexyplop *)

(* begin module binomial.themain *)
procedure themain(var xyin, xyplop, binomialp: text);
(* the main procedure of the binomial program *)
var
   cumulative: real; (* cumulative probabilities *)
   black, white: integer; (* number of black and white balls found *)
   p: real; (* the probability of black balls *)
   q: real; (* the probability of white balls *)
   maxwhite: integer; (* maximum number of blacks to plot *)
   mean: real; (* the mean of the distribution *)
   pbin: real; (* the calculated probability *)
   sigma: real; (* the standard deviation of the distribution *)
   n: integer; (* the total of blacks and whites, number of samples *)
begin
   writeln(output,'binomial ',version:4:2);

   reset(binomialp);
   readln(binomialp,n);
   readln(binomialp,p);
   readln(binomialp,maxwhite);

writeln (output,'n = ',n:1);
writeln (output,'p = ',p:10:8);
writeln (output,'maxwhite = ',maxwhite:1);

   if (p < 0) or (p > 1) then begin
      writeln(output,'probability must be between 0 and 1.');
      halt
   end;

   if black < 0 then begin
      writeln(output,'there cannot be negative black balls');
      halt
   end;

   if white < 0 then begin
      writeln(output,'there cannot be negative white balls');
      halt
   end;

{
old 
   n := black + white;
}

   if n < 0 then begin
      writeln(output,'there must be a positive number of balls');
      halt
   end;

   q := 1-p;
   mean := n*p;
   sigma := sqrt(n*p*q);
   if maxwhite > n then maxwhite := n;

   rewrite(xyplop);
   makexyplop(xyplop,maxwhite);

   rewrite(xyin);
   writeln(xyin,'* binomial ',version:4:2);
   writeln(xyin,'* n=',n:1);
   writeln(xyin,'* probability of white = ',p:5:3);
   writeln(xyin,'* highest white value  = ',maxwhite:4);
   writeln(xyin,'* -----');
   writeln(xyin,'* mean  = ',mean:4:2);
   writeln(xyin,'* sigma = ',sigma:4:2);
   writeln(xyin,'* column 1: number of white balls');
   writeln(xyin,'* column 2: probability of withdrawing that number white balls in');
   writeln(xyin,'*           ',n:1,' black and white balls total');
   writeln(xyin,'*           given a probability of ',p:5:3,
                ' white balls in the population');
   writeln(xyin,'* column 3: cumulative sum of column 2');

   writeln(output,'                  mean = ',mean:5:2);
   writeln(output,'                 sigma = ',sigma:5:2);
   writeln(output,'sharpness = mean/sigma = ',mean/sigma:5:2);

   cumulative := 0.0;
   for white := 0 to maxwhite do begin
      black := n - white;
      pbin := pbinomial(black,white,p);
      cumulative := pbin + cumulative; 
      writeln(xyin,white:8,' ',pbin:15:13,' ', cumulative:15:13)
   end
end;
(* end module binomial.themain *)

begin (* binomial *)
   themain(xyin, xyplop, binomialp);
1: end.