/* Output from p2c 1.21alpha-07.Dec.93, the Pascal-to-C translator */ /* From input file "binomial.p" */ #include /* 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 */ /* end of program */ /* begin module version */ #define 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 */ Static _TEXT xyin; /* the input file to xyplop */ Static _TEXT xyplop; /* the control file for xyplop */ Static _TEXT binomialp; /* control parameters for binomial */ Static jmp_buf _JL1; /* begin module halt */ Static Void 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. */ printf(" program halt.\n"); longjmp(_JL1, 1); } Local double combination(black, white) long black, white; { /* calculate the combinations for a number of black and white balls. give the result as the log of the power */ long index; /* counting index */ double sum = 0.0; /* the sum of the logs of the factorial components */ /* 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; index <= black + white; index++) sum += log((double)index); /* now remove the factorial of white */ for (index = 1; index <= white; index++) sum -= log((double)index); return sum; } /* combination */ Local double power(black, white, fraction) long black, white; double fraction; { /* calculate the power portion of the binomial distribution. give the result as the log of the power */ if (fraction >= 1.0) return 0.0; else return (black * log(fraction) + white * log(1 - fraction)); } /* power */ /* end module halt version = 'delmod 6.16 84 mar 12 tds/gds'; */ /* begin module pbinomial */ Static double pbinomial(black, white, fraction) long black, white; double fraction; { /* 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. */ double sum; /* log of the solution */ long total; /* sum of black and white */ 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) return 0.0; else return exp(sum); } /* end module pbinomial */ /* begin module tpbinom */ Static Void tpbinom(black, white, fraction) long black, white; double fraction; { /* test the pbinomial procedure for the given fraction */ long index; /* counter for the running sum */ double sum = 0.0; /* the running sum of the probabilities */ long total; /* sum of black and white */ printf("testing pbinomial procedure\n"); total = black + white; for (index = 0; index <= total; index++) sum += pbinomial(index, total - index, fraction); printf("sum of p(b,w|fraction) for %ld balls is %10.8f\n", total, sum); } /* end module tpbinom */ /* begin module binomial.makexyplop */ Static Void makexyplop(x, maxwhite) _TEXT *x; long maxwhite; { /* write the xyplop values out */ long xintervals; /* for computing intervals to show */ if (maxwhite > 20) xintervals = 10; else xintervals = maxwhite; fprintf(x->f, "3 6 zerox zeroy\n"); fprintf(x->f, "zx 0 %ld zx min max \n", maxwhite); fprintf(x->f, "zy 0 1 zy min max \n"); fprintf(x->f, "%ld 10 1 1 xinterval yinterval \n", xintervals); fprintf(x->f, "6 6 xwidth ywidth \n"); fprintf(x->f, "0 3 xdecimal ydecimal \n"); fprintf(x->f, "15 15 xsize ysize \n"); fprintf(x->f, "number of white\n"); /* label x */ fprintf(x->f, "p(number of white)\n"); /* label y */ /*writeln(x,'p(b=',black:1,': w=',white:1, ' | number b)'); old label */ fprintf(x->f, "a no cross hairs\n"); fprintf(x->f, "n 2 no log\n"); fprintf(x->f, "n 2 no log\n"); fprintf(x->f, " *-\n"); fprintf(x->f, "1 2 xcolumn ycolumn\n"); fprintf(x->f, "0 symbol column\n"); fprintf(x->f, "0 0 xscolumn yscolumn\n"); fprintf(x->f, "0 0 0 hue saturation brightness columns for color manipulation\n"); fprintf(x->f, " *-\n"); fprintf(x->f, "c symbol-to-plot\n"); fprintf(x->f, "o symbol flag\n"); fprintf(x->f, "0.250 symbol x size\n"); fprintf(x->f, "0.250 symbol y size\n"); fprintf(x->f, "cl 0.05 connected line \n"); fprintf(x->f, "n 0.05 no regression-line \n"); fprintf(x->f, " *-\n"); fprintf(x->f, ".\n"); fprintf(x->f, " *-\n"); fprintf(x->f, " *-\n"); fprintf(x->f, "p 1.50 0.50 1.50 1.00 edgecontrol (p=page), edgeleft, edgeright, edgelow, edgehigh in cm\n"); fprintf(x->f, "9.01 version of xyplo that this parameter file is designed for.\n"); } /* end module binomial.makexyplop */ /* begin module binomial.themain */ Static Void themain(xyin, xyplop, binomialp) _TEXT *xyin, *xyplop, *binomialp; { /* the main procedure of the binomial program */ double cumulative = 0.0; /* cumulative probabilities */ long black, white; /* number of black and white balls found */ double p; /* the probability of black balls */ double q; /* the probability of white balls */ long maxwhite; /* maximum number of blacks to plot */ double mean; /* the mean of the distribution */ double pbin; /* the calculated probability */ double sigma; /* the standard deviation of the distribution */ long n; /* the total of blacks and whites, number of samples */ printf("binomial %4.2f\n", version); if (*binomialp->name != '\0') { if (binomialp->f != NULL) binomialp->f = freopen(binomialp->name, "r", binomialp->f); else binomialp->f = fopen(binomialp->name, "r"); } else rewind(binomialp->f); if (binomialp->f == NULL) _EscIO2(FileNotFound, binomialp->name); RESETBUF(binomialp->f, Char); fscanf(binomialp->f, "%ld%*[^\n]", &n); getc(binomialp->f); fscanf(binomialp->f, "%lg%*[^\n]", &p); getc(binomialp->f); fscanf(binomialp->f, "%ld%*[^\n]", &maxwhite); getc(binomialp->f); printf("n = %ld\n", n); printf("p = %10.8f\n", p); printf("maxwhite = %ld\n", maxwhite); if ((unsigned)p > 1) { printf("probability must be between 0 and 1.\n"); halt(); } if (black < 0) { printf("there cannot be negative black balls\n"); halt(); } if (white < 0) { printf("there cannot be negative white balls\n"); halt(); } /* old n := black + white; */ if (n < 0) { printf("there must be a positive number of balls\n"); halt(); } q = 1 - p; mean = n * p; sigma = sqrt(n * p * q); if (maxwhite > n) maxwhite = n; if (*xyplop->name != '\0') { if (xyplop->f != NULL) xyplop->f = freopen(xyplop->name, "w", xyplop->f); else xyplop->f = fopen(xyplop->name, "w"); } else { if (xyplop->f != NULL) rewind(xyplop->f); else xyplop->f = tmpfile(); } if (xyplop->f == NULL) _EscIO2(FileNotFound, xyplop->name); SETUPBUF(xyplop->f, Char); makexyplop(xyplop, maxwhite); if (*xyin->name != '\0') { if (xyin->f != NULL) xyin->f = freopen(xyin->name, "w", xyin->f); else xyin->f = fopen(xyin->name, "w"); } else { if (xyin->f != NULL) rewind(xyin->f); else xyin->f = tmpfile(); } if (xyin->f == NULL) _EscIO2(FileNotFound, xyin->name); SETUPBUF(xyin->f, Char); fprintf(xyin->f, "* binomial %4.2f\n", version); fprintf(xyin->f, "* n=%ld\n", n); fprintf(xyin->f, "* probability of white = %5.3f\n", p); fprintf(xyin->f, "* highest white value = %4ld\n", maxwhite); fprintf(xyin->f, "* -----\n"); fprintf(xyin->f, "* mean = %4.2f\n", mean); fprintf(xyin->f, "* sigma = %4.2f\n", sigma); fprintf(xyin->f, "* column 1: number of white balls\n"); fprintf(xyin->f, "* column 2: probability of withdrawing that number white balls in\n"); fprintf(xyin->f, "* %ld black and white balls total\n", n); fprintf(xyin->f, "* given a probability of %5.3f white balls in the population\n", p); fprintf(xyin->f, "* column 3: cumulative sum of column 2\n"); printf(" mean = %5.2f\n", mean); printf(" sigma = %5.2f\n", sigma); printf("sharpness = mean/sigma = %5.2f\n", mean / sigma); for (white = 0; white <= maxwhite; white++) { black = n - white; pbin = pbinomial(black, white, p); cumulative = pbin + cumulative; fprintf(xyin->f, "%8ld %15.13f %15.13f\n", white, pbin, cumulative); } } /* end module binomial.themain */ main(argc, argv) int argc; Char *argv[]; { PASCAL_MAIN(argc, argv); if (setjmp(_JL1)) goto _L1; binomialp.f = NULL; strcpy(binomialp.name, "binomialp"); xyplop.f = NULL; strcpy(xyplop.name, "xyplop"); xyin.f = NULL; strcpy(xyin.name, "xyin"); themain(&xyin, &xyplop, &binomialp); _L1: if (xyin.f != NULL) fclose(xyin.f); if (xyplop.f != NULL) fclose(xyplop.f); if (binomialp.f != NULL) fclose(binomialp.f); exit(EXIT_SUCCESS); } /* End. */