/* Output from p2c 1.21alpha-07.Dec.93, the Pascal-to-C translator */ /* From input file "binom.p" */ #include /* produce the binom probabilities for a found black to white ratio Thomas Schneider */ /* end of program */ /* begin module version */ #define 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 */ Static _TEXT xyin; /* the input file to xyplop */ Static _TEXT xyplop; /* the control file for xyplop */ Static _TEXT binomp; /* control parameters for binom */ 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, P) long black, white; double P; { /* calculate the power portion of the binomial distribution. give the result as the log of the power */ if (P >= 1.0) return 0.0; else return (black * log(P) + white * log(1 - P)); } /* power */ /* end module halt version = 'delmod 6.16 84 mar 12 tds/gds'; */ /* begin module pbinomial */ Static double pbinomial(black, white, P) long black, white; double P; { /* 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. */ double sum; /* log of the solution */ long total; /* sum of black and white */ 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) return 0.0; else return exp(sum); } /* end module pbinomial */ /* begin module tpbinom */ Static Void tpbinom(black, white, P) long black, white; double P; { /* test the pbinomial procedure for the given P */ 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, P); printf("sum of p(b,w|P) for %ld balls is %10.8f\n", total, sum); } /* end module tpbinom */ /* begin module binom.makexyplop */ Static Void makexyplop(x, maxblack) _TEXT *x; long maxblack; { /* write the xyplop values out */ fprintf(x->f, "1 2 zerox zeroy\n"); fprintf(x->f, "x 0 %ld zx min max \n", maxblack); fprintf(x->f, "y 0 1 zy min max \n"); fprintf(x->f, "20 10 1 1 xinterval yinterval \n"); fprintf(x->f, "4 5 xwidth ywidth \n"); fprintf(x->f, "0 3 xdecimal ydecimal \n"); fprintf(x->f, "6 6 xsize ysize \n"); fprintf(x->f, "number of blacks\n"); /* label x */ fprintf(x->f, "p(number of blacks)\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 no log\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, "0 0 0 hue saturation brightness columns for color manipulation\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 binom.makexyplop */ /* begin module binom.themain */ Static Void themain(xyin, xyplop, binomp) _TEXT *xyin, *xyplop, *binomp; { /* the main procedure of the binom program */ long black, white; /* number of black and white balls found */ double P; /* the probability of black balls */ double Q; /* the probability of white balls */ double mean; /* the mean of the distribution */ double pbin; /* the calculated probability */ double sigma; /* the standard deviation of the distribution */ long total; /* the total of blacks and whites */ printf("binom %4.2f\n", version); if (*binomp->name != '\0') { if (binomp->f != NULL) binomp->f = freopen(binomp->name, "r", binomp->f); else binomp->f = fopen(binomp->name, "r"); } else rewind(binomp->f); if (binomp->f == NULL) _EscIO2(FileNotFound, binomp->name); RESETBUF(binomp->f, Char); fscanf(binomp->f, "%ld%ld%*[^\n]", &black, &white); getc(binomp->f); fscanf(binomp->f, "%lg%*[^\n]", &P); getc(binomp->f); if (P < 0) { printf("probability P must be positive\n"); halt(); } if (P > 1) { printf("probability P must be less than 1.0\n"); halt(); } total = black + white; Q = 1 - P; mean = total * P; sigma = sqrt(total * P * Q); /* rewrite(xyplop); */ 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, "* binom %4.2f\n", version); fprintf(xyin->f, "* black=%ld: white=%ld\n", black, white); fprintf(xyin->f, "* probability of black =%5.3f\n", P); fprintf(xyin->f, "* mean = %4.2f\n", mean); fprintf(xyin->f, "* sigma = %4.2f\n", sigma); fprintf(xyin->f, "* Zblack = %4.2f\n", (black - mean) / sigma); /* 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); fprintf(xyin->f, "probability of this black/white combination: %10.8f\n", pbin); } /* end module binom.themain */ main(argc, argv) int argc; Char *argv[]; { PASCAL_MAIN(argc, argv); if (setjmp(_JL1)) goto _L1; binomp.f = NULL; strcpy(binomp.name, "binomp"); xyplop.f = NULL; strcpy(xyplop.name, "xyplop"); xyin.f = NULL; strcpy(xyin.name, "xyin"); themain(&xyin, &xyplop, &binomp); _L1: if (xyin.f != NULL) fclose(xyin.f); if (xyplop.f != NULL) fclose(xyplop.f); if (binomp.f != NULL) fclose(binomp.f); exit(EXIT_SUCCESS); } /* End. */