/* Output from p2c 1.21alpha-07.Dec.93, the Pascal-to-C translator */ /* From input file "binplo.p" */ #include /* produce the binomial probabilities for a found black to white ratio Thomas Dana Schneider 1987 */ /* end of program */ /* begin module version */ #define 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 */ Static _TEXT xyin; /* the input file to xyplop */ Static _TEXT xyplop; /* the control file for xyplop */ Static _TEXT binplop; /* control parameters for binplo */ 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 */ 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. */ long total; /* sum of black and white */ total = black + white; return exp(combination(black, white) + power(black, white, fraction)); } /* 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 binplo.makexyplop */ Static Void makexyplop(x, black, white) _TEXT *x; long black, white; { /* write the xyplop values out */ fprintf(x->f, "x 0 1 zx min max \n"); fprintf(x->f, "y 0 %10.8f zy min max \n", pbinomial(black, white, (double)black / (black + white))); fprintf(x->f, "20 20 xinterval yinterval \n"); fprintf(x->f, "4 5 xwidth ywidth \n"); fprintf(x->f, "2 3 xdecimal ydecimal \n"); fprintf(x->f, "5.5 4 xsize ysize \n"); fprintf(x->f, "fraction of black\n"); /* label x */ fprintf(x->f, "p(b=%ld: w=%ld | fraction b)\n", black, white); /* label y */ fprintf(x->f, "n no cross hairs\n"); fprintf(x->f, "n 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, " --\n"); fprintf(x->f, "+ symbol-to-plot\n"); fprintf(x->f, "+ symbol flag\n"); fprintf(x->f, "0.01 symbol x size\n"); fprintf(x->f, "0.01 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"); } /* end module binplo.makexyplop */ /* begin module binplo.themain */ Static Void themain(xyin, xyplop, binplop) _TEXT *xyin, *xyplop, *binplop; { /* the main procedure of the binplo program */ long black, white; /* number of black and white balls found */ double fraction = 0.0; /* the true fraction of black and white balls */ long index; /* index to the plotting of the points */ double interval; /* 1/index */ long points; /* number of points to plot */ double p; /* the calculated probability */ printf("binplo %4.2f\n", version); if (*binplop->name != '\0') { if (binplop->f != NULL) binplop->f = freopen(binplop->name, "r", binplop->f); else binplop->f = fopen(binplop->name, "r"); } else rewind(binplop->f); if (binplop->f == NULL) _EscIO2(FileNotFound, binplop->name); RESETBUF(binplop->f, Char); fscanf(binplop->f, "%ld%ld%*[^\n]", &black, &white); getc(binplop->f); fscanf(binplop->f, "%ld%*[^\n]", &points); getc(binplop->f); 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, black, white); 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, "* binplo %4.2f\n", version); fprintf(xyin->f, "* %ld intervals plotted\n", points); fprintf(xyin->f, "* black=%ld: white=%ld\n", black, white); if (points <= 0) { if (points == 0) tpbinom(black, white, (double)black / (black + white)); else { printf(" you cannot have a negative number of points\n"); halt(); } return; } interval = 1.0 / points; for (index = 1; index <= points; index++) { fraction += interval; p = pbinomial(black, white, fraction); if ((long)(p * 1e8) > 0) fprintf(xyin->f, "%10.8f %10.8f\n", fraction, p); } } /* end module binplo.themain */ main(argc, argv) int argc; Char *argv[]; { PASCAL_MAIN(argc, argv); if (setjmp(_JL1)) goto _L1; binplop.f = NULL; strcpy(binplop.name, "binplop"); xyplop.f = NULL; strcpy(xyplop.name, "xyplop"); xyin.f = NULL; strcpy(xyin.name, "xyin"); themain(&xyin, &xyplop, &binplop); _L1: if (xyin.f != NULL) fclose(xyin.f); if (xyplop.f != NULL) fclose(xyplop.f); if (binplop.f != NULL) fclose(binplop.f); exit(EXIT_SUCCESS); } /* End. */