/* Output from p2c 1.21alpha-07.Dec.93, the Pascal-to-C translator */ /* From input file "multtest.p" */ #include /* multtest: multiple Student's t-test Dr. Thomas D. Schneider National Institutes of Health National Cancer Institute Center for Cancer Research Nanobiology Program Molecular Information Theory Group Frederick, Maryland 21702-1201 toms@ncifcrf.gov permanent email: toms@alum.mit.edu (use only if first address fails) http://www.ccrnp.ncifcrf.gov/~toms/ */ /* end of program */ /* begin module version */ #define version 1.21 /* of multtest.p 2007 Mar 21 2007 Mar 21, 1.21: upgrade simpson 2001 Jun 14, 1.20: documentation upgrade 1.19, 2000 April 19: working version origin from version 1.07 of ttest.p 1996 September 30 origin of ttest was probably 1994 January 15 */ /* end module version */ /* begin module describe.multtest */ /* name multtest: multiple Student's t-test synopsis multtest(distributions: in, multtestp: in, list: out, xyin: out, key: out, output: out) files distributions: parameters to control the program: A set of 3 numbers on one line defines a distribution: N (integer) mean (real) standard deviation (real) This file contains a series of such distributions, one per line. fudgefactor: If a line begins with 'f' then all standard deviations from that point on (until the next 'f') will be multiplied by the number following on the same line as f. This allows one to see the effects if the standard deviation does not account for all uncertainties. Lines that begin with '*' are comments. multtestp: a single character OR pairs of numbers, one per line, that represent the distributions to be compared by a ttest. The distributions listed are the only ones output to the list file. controlchar: If a line begins with a character then: l: compute the natural log of the mean and adjust the standard deviation accordingly (divide by the mean). c: set the color transform. The next character is the colortype: the colors produced in xyin are adjusted: n: normal colors Hue is set to the probability. i: interval colors. The 'i' must be followed by the number of intervals to truncate the colors to. Hue is set to the truncated probability. For example, if the number of intervals is 2, then the colors will be 0.0 to 0.5 and 0.5 to 1. s: standard interval colors. The colors will cut off at the points 0.50, 0.90 to 0.95 and 0.99. list: Input values, calculated Student's t value and probabilities that the numbers are the same, along with the Z score for the energy deviation from zero, and the probability that this deltadeltaG is zero. The pairs to be listed are defined in the multtestp file. xyin: input to xyplo program: column 1: (integer) first distribution number column 2: (integer) second distribution number column 3: (integer) hue fraction column 4: (integer) saturation fraction = 1 column 5: (integer) brightness fraction = 1 The hue is computed from the probability according to the colortype defined in multtestp. key: input to xyplo program, same colomuns but generates the key for the probabilities. output: messages to the user description This program performs T test computations for pairs of distributions and converts them to probabilities that the means are the same. This is a two-tailed Student's t test. examples For a multtestp file: 1 2 with a distributions file: 3 118 36 3 205 26 The resulting list file is: multtest 1.02 distribution 1 | distribution 2 number 3 | 3 mean 118.00000 | 205.00000 standard dev. 36.00000 | 26.00000 sigma-D = 25.63851 degrees of freedom = 4 t = 3.39333 p(1=2) = 0.02745 This is the probability that they are the same. documentation @book{Press1989, author = "W. H. Press and B. P. Flannery and S. A. Teukolsky and W. T. Vetterling", title = "Numerical Recipies in Pascal. The Art of Scientific Computing", publisher = "Cambridge University Press", address = "Cambridge", year = "1989"} Given a t value from a Student's t test, and the degrees of freedom, df, return the probability for a two tailed test. The code for doing this was originally in java script, from: Richard Lowry Department of Psychology Vassar College Poughkeepsie, NY 12604-0396 USA office: (914)437-7381 fax: (914)437-7538 lowry@vassar.edu http://faculty.vassar.edu/~lowry/VassarStats.html The original functional html containing this code is also given. It was translated to Pascal by Tom Schneider. see also {The source program:} ttest.p {The xyin file can be graphed using:} xyplo.p {Results can be checked at:} http://www.statsol.com/tools/stattools/ttestindependenttool.html http://faculty.vassar.edu/~lowry/VassarStats.html {The method of conversion to a log scale comes:} calc.p author Thomas Dana Schneider bugs technical notes */ /* end module describe.multtest */ /* begin module calc.const */ #define defrwidth 9 /* default width of numbers */ #define defrdecimal 2 /* default decimal places of numbers */ #define zbound 38 /* largest z that the g function can tolerate without bombing The value zbound = 38 works for a Sun Sparc workstation. */ #define pmstring1 " \\pm " /* latex format */ #define pmstring2 " +/- " /* regular format */ /* end module calc.const */ /* from calc.p */ typedef struct range { /* representation of x +/- error */ double estimate; /* the central value x */ double error; /* the most probable values above and below x */ } range; Static _TEXT distributions, multtestp, list, xyin, key; /* files used by this program */ /* calc variables */ Static boolean showingerrors; /* true if error ranges are shown */ /* not used: latex: boolean; (* if true show "\pm" otherwise use "+/-" for range. *) */ Static Char pmstring[5]; /* string to show depending on latex variable */ Static boolean fromto; Static jmp_buf _JL1; /* true = [5 to 11] display, false = 8 +/- 3 display */ /* 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 variables for writerange: */ struct LOC_writerange { _TEXT *fout; range r; long width, decimal; boolean exponential, zerodecimal; /* true if there are to be no decimal places */ } ; Local Void makeplusminus(LINK) struct LOC_writerange *LINK; { if (LINK->exponential) { fprintf(LINK->fout->f, "% .*E%.5s% .*E", P_max((int)LINK->width - 7, 1), LINK->r.estimate, pmstring, P_max((int)LINK->width - 7, 1), LINK->r.error); return; } if (LINK->zerodecimal) fprintf(LINK->fout->f, "%*ld%.5s%*ld", (int)LINK->width, (long)floor(LINK->r.estimate + 0.5), pmstring, (int)LINK->width, (long)floor(LINK->r.error + 0.5)); else { fprintf(LINK->fout->f, "%*.*f%.5s%*.*f", (int)LINK->width, (int)LINK->decimal, LINK->r.estimate, pmstring, (int)LINK->width, (int)LINK->decimal, LINK->r.error); /* then write(fout, r.estimate:width ,' +/- ', r.error:width ) else if zerodecimal then write(fout, round(r.estimate):width,' +/- ', round(r.error):width) else write(fout, r.estimate:width:decimal,' +/- ', r.error:width:decimal) */ } } Local Void makefromto(LINK) struct LOC_writerange *LINK; { putc('[', LINK->fout->f); if (LINK->exponential) fprintf(LINK->fout->f, "% .*E to % .*E", P_max((int)LINK->width - 7, 1), LINK->r.estimate - LINK->r.error, P_max((int)LINK->width - 7, 1), LINK->r.estimate + LINK->r.error); else if (LINK->zerodecimal) fprintf(LINK->fout->f, "%*ld to %*ld", (int)LINK->width, (long)floor(LINK->r.estimate - LINK->r.error + 0.5), (int)LINK->width, (long)floor(LINK->r.estimate + LINK->r.error + 0.5)); else fprintf(LINK->fout->f, "%*.*f to %*.*f", (int)LINK->width, (int)LINK->decimal, LINK->r.estimate - LINK->r.error, (int)LINK->width, (int)LINK->decimal, LINK->r.estimate + LINK->r.error); putc(']', LINK->fout->f); } /* end module halt version = 'matmod 2.07 2007mar21 tds/gds'; */ /* begin module writerange */ Static Void writerange(fout_, r_, width_, decimal_, exponential_) _TEXT *fout_; range r_; long width_, decimal_; boolean exponential_; { /* write the range to fout, using the given width and decimal places for the numbers. If exponential is true, then use exponential form */ struct LOC_writerange V; V.fout = fout_; V.r = r_; V.width = width_; V.decimal = decimal_; V.exponential = exponential_; V.zerodecimal = (V.decimal <= 0); if (showingerrors) { if (fromto) makefromto(&V); else makeplusminus(&V); return; } if (V.exponential) { fprintf(V.fout->f, "% .*E", P_max((int)V.width - 7, 1), V.r.estimate); return; } if (V.zerodecimal) fprintf(V.fout->f, "%*ld", (int)V.width, (long)floor(V.r.estimate + 0.5)); else fprintf(V.fout->f, "%*.*f", (int)V.width, (int)V.decimal, V.r.estimate); } /* writerange */ /* end module writerange */ /* begin module timesrange */ Static Void timesrange(a, b, atimesb) range a, b, *atimesb; { /* multiply a by b and put result in atimesb */ double TEMP, TEMP1; atimesb->estimate = a.estimate * b.estimate; if (a.estimate != 0 && b.estimate != 0) { TEMP = a.error / a.estimate; TEMP1 = b.error / b.estimate; atimesb->error = fabs(atimesb->estimate) * sqrt(TEMP * TEMP + TEMP1 * TEMP1); return; } printf("WARNING: one of the estimates to be multiplied is zero\n"); printf("so the error will be set to zero\n"); atimesb->error = 0.0; } /* end module timesrange */ /* begin module dividerange */ Static Void dividerange(a, b, aoverb) range a, b, *aoverb; { /* divide a by b and put result in aoverb */ double TEMP, TEMP1; aoverb->estimate = a.estimate / b.estimate; if (a.estimate != 0) { TEMP = a.error / a.estimate; TEMP1 = b.error / b.estimate; aoverb->error = fabs(aoverb->estimate) * sqrt(TEMP * TEMP + TEMP1 * TEMP1); return; } printf("WARNING: one of the estimates to be divided is zero\n"); printf("so the error will be set to zero\n"); aoverb->error = 0.0; } /* end module dividerange */ /* begin module logrange */ Static Void logrange(a, base, loga) range a; double base; range *loga; { /* take log base of a as loga */ loga->estimate = log(a.estimate) / log(base); loga->error = a.error / a.estimate; } /* end module logrange */ /* begin module simpson */ Static Void simpson(upper, answer, tol) double upper, *answer, *tol; { /* Perform a numerical integration of the normal distribution using Simpson's rule. The variable "lower" is the lower bound of the integration. The variable "upper" is the the upper bound of the integration. "answer" is the value of the integration, and "tol" is the tolerance of the integration procedure, and thus the error of the calculation. Written by Mark Shaner, 1992. Taken from Miller, Alan R. PASCAL PROGRAMS FOR SCIENTISTS AND ENGINEERS. SYBEX Inc., 1981. pgs 260-291 2007 Mar 21: Tom Schneider: Using a tolerance of 1/maxint caused an infinite loop on a 64 bit machine. Presumably it was attempting to compute the value to 20 or so decimal places. To avoid this, the tolerance is set to a more reasonable value. */ /* a counter */ long i; double x, x1, x2; /* the independent variables for calculating the value of functions */ double pi = 4.0 * atan(1.0); /* the value of pi */ double val; /* the value of 1/sqrt(2*pi), it is used in calculating the value of the gaussian for any x value, and is defined as a variable in order to speed calculations. */ double deltax; /* the distance between every x value and the subsequent one */ double evensum = 0.0; /* the sum of the area under each of the even numbered parabolas */ double oddsum; /* the sum of the area under each of the odd numbered parabolas */ double endsum; /* the sum of the area under the first and last parabolas */ double endcor; /* the value of the end correction, it is determined using dgauss */ double sum1; /* a place to store the previous sum so that it can be compared with the subsequent sum to determine if the tolerance level has been reached */ long pieces = 2; /* the number of parabolas under the curve */ double lower = 0.0; /* the lower bound of the integration */ double sum; /* the value of the area under the curve */ val = 1 / sqrt(2 * pi); /* activate both files to be used */ upper = fabs(upper); /* tol := 1/maxint; */ *tol = 1.0e-6; deltax = (upper - lower) / pieces; x = lower + deltax; oddsum = val * exp(-0.5 * x * x); x1 = lower; x2 = upper; endsum = val * exp(-0.5 * x1 * x1) + val * exp(-0.5 * x2 * x2); endcor = x2 * val * exp(-0.5 * x2 * x2) - x1 * val * exp(-0.5 * x1 * x1); sum = (endsum + 4.0 * oddsum) * deltax / 3.0; do { pieces *= 2; sum1 = sum; deltax = (upper - lower) / pieces; evensum += oddsum; oddsum = 0.0; for (i = 1; i <= pieces / 2; i++) { x = lower + deltax * (2.0 * i - 1.0); oddsum += val * exp(-0.5 * x * x); } sum = (7.0 * endsum + 14.0 * evensum + 16.0 * oddsum + endcor * deltax) * deltax / 15.0; } while (fabs(sum - sum1) > fabs(*tol * sum)); *answer = 0.5 - sum; if (*answer < 0.0) /* safety for roundoff errors */ *answer = 0.0; /* writeln (upper:7:5,' ',answer:7:5,'+/-',tol) */ } /* Local variables for ttestprobability: */ struct LOC_ttestprobability { double pj2; /* pi/2 */ } ; Local double zip(q, i, j, b, LINK) double q, i, j, b; struct LOC_ttestprobability *LINK; { double k; double zz = 1.0; double z; z = zz; k = i; while (k <= j) { zz = zz * q * k / (k - b); z += zz; k += 2; } return z; } Local double buzz(t, n, LINK) double t, n; struct LOC_ttestprobability *LINK; { double rt, fk, ek, dk; t = fabs(t); rt = t / sqrt(n); fk = atan(rt); if (n == 1) return (1 - fk / LINK->pj2); else { ek = sin(fk); dk = cos(fk); if ((((long)floor(n + 0.5)) & 1) == 1) return (1 - (fk + ek * dk * zip(dk * dk, 2.0, n - 3, -1.0, LINK)) / LINK->pj2); else return (1 - ek * zip(dk * dk, 1.0, n - 3, -1.0, LINK)); } } /* end module simpson version = 'matmod 2.07 2007mar21 tds/gds'; */ /******************************************************************************/ /* begin module ttestprobability */ Static double ttestprobability(t, df) double t, df; { /* Given a t value from a Student's t test, and the degrees of freedom, df, return the probability for a two tailed test. The code was originally in java script, from: Richard Lowry Department of Psychology Vassar College Poughkeepsie, NY 12604-0396 USA office: (914)437-7381 fax: (914)437-7538 lowry@vassar.edu http://faculty.vassar.edu/~lowry/VassarStats.html The original functional html containing this code is given below the Pascal. It was translated to Pascal by: Dr. Thomas D. Schneider National Cancer Institute Laboratory of Experimental and Computational Biology Frederick, Maryland 21702-1201 toms@ncifcrf.gov permanent email: toms@alum.mit.edu http://www.lecb.ncifcrf.gov/~toms/ It is tested carefully in: testttestprobability. */ struct LOC_ttestprobability V; double pi = 4.0 * atan(1.0); /* ration of circle circumfrence to diameter */ double pj4; /* pi/2 */ double pi2; /* 2*pi */ double e = exp(1.0); /* e */ double dgr; /* 180/pi */ /* NOT USED function Abuzz(p,n: real): real; var v, dv, t: real; begin v := 0.5; dv := 0.5; t := 0; while (dv>1e-6) do begin t := 1/v-1; dv := dv/2; if(buzz(t,n)>p) then v := v-dv else v := v+dv end; Abuzz := t end; */ V.pj2 = pi / 2; pj4 = pi / 4; pi2 = 2 * pi; dgr = 180 / pi; return (1.0 - buzz(t, df, &V) / 2.0); } /* ttestprobability */ #define wid 6 /* width of output */ #define dec 3 /* decimals of output */ Local Void try_(t, df) double t, df; { /* compute and show a test t and df value t: t value, df: degrees of freedom */ printf(" t = %*.*f df = %*ld P = %*.*f\n", wid, dec, t, wid, (long)floor(df + 0.5), wid, dec, ttestprobability(t, df)); } /* try */ /* The material below this point is a functioning javascript program that can be used as a web page. ******************************************************************************** t to P

version = 2.00 of t2p.html 2000 April 18
  t =
df =
  P (one-tailed)
  P (two-tailed)
p * (upperhue - lowerhue) + lowerhue; p2 = p1 / 0.85; /* color adjusted */ switch (LINK->colortype) { case 'n': /* normal */ LINK->hue = p2; break; case 'i': /* intervals */ LINK->hue = (double)((long)(LINK->intervals * p2)) / LINK->intervals; break; case 's': /* standard intervals */ if (p2 < 0.50) LINK->hue = color1; else if (p2 < 0.90) LINK->hue = color2; else if (p2 < 0.95) LINK->hue = color3; else if (p2 < 0.99) LINK->hue = color4; else { LINK->hue = color5; /* OLD STUFF: if p < 0.01 then hue := color1 else if p < 0.05 then hue := color2 else if p < 0.10 then hue := color3 else hue := color4 if p < 0.01 then hue := color1 else if p < 0.05 then hue := color2 else if p < 0.10 then hue := color3 else hue := color4 */ /* if p < 0.05 then hue := 0.0 else if p < 0.10 then hue := 0.5 else hue := 1.0 */ } break; } /* saturation := 0.5; brightness := 1.0; brightness := (1.0-p*p*p); saturation := p; if p < 0.1 then saturation := p else saturation := 0.5; if p > 0 then saturation := 1/p else saturation := 0.0; saturation := 0.8; if p > 0 then saturation := 1/(p*p) else saturation := 1.0; if p < 0.1 then if p > 0 then saturation := p else saturation := 0.0 else saturation := p; if p < 0.1 then saturation := 0.0 else saturation := p; 2000 Apr 4 Kd color scheme: saturation := 10*p; if saturation > 1.0 then saturation := 1.0; (* set saturation to 1 if p > 0.02, otherwise ramp down to zero. This sets the lowest color to white. *) saturation := 500000*p*p*p; (* 1% cutoff method:*) if p > 0.010 then saturation := 1.0 else saturation := 0.1 (* show probability of inequality *) if p < 0.02 then begin hue := p*50; saturation := 1.0; brightness := 1.0; end; else begin hue := 1.0; saturation := 0.0; brightness := 1.0; end; (* 0.02 cutoff, used with energy fig 2000 Apr 4 *) saturation := 50*p; if saturation > 1.0 then saturation := 1.0; (* 0.02 cutoff *) saturation := 50*p*p; if p > 0.05 then saturation := 1.0; if saturation > 1.0 then saturation := 1.0; */ /* 95 % confidence limits */ if (LINK->p > 1 - 0.95) { /* do colors */ LINK->saturation = 1.0; LINK->brightness = 1.0; LINK->symbol = 'c'; return; } /* do hollowsquare */ LINK->saturation = 0.0; LINK->brightness = 0.0; LINK->symbol = 'b'; } #undef lowerhue #undef upperhue Local Void showcolor(xyin, d1, d2, symbol, LINK) _TEXT *xyin; double d1, d2; Char symbol; struct LOC_themain *LINK; { /* show the color to xyin at coordinate (d1, d2).*/ fprintf(xyin->f, "%*.*f %*.*f %*.*f %*.*f %*.*f", wid, dec, d1 + colorshift, wid, dec, d2 + colorshift, wid, dec, LINK->hue, wid, dec, LINK->saturation, wid, dec, LINK->brightness); if (symbol == 'c') fprintf(xyin->f, " %c", symbol); else fprintf(xyin->f, " hollowsquare"); putc('\n', xyin->f); /* writeln(xyin, (d1+colorshift):wid:dec, ' ',(d2+colorshift):wid:dec, ' ',hue:wid:dec, ' ',saturation:wid:dec, ' ',brightness:wid:dec); writeln(xyin, d1:1, ' ',d2:1, ' ',hue:wid:dec, ' ',saturation:wid:dec, ' ',brightness:wid:dec); */ } Local Void compute(LINK) struct LOC_themain *LINK; { /* compute the t and p values using the global variables */ double TEMP, TEMP1; LINK->df = LINK->n[LINK->d1-1] + LINK->n[LINK->d2-1] - 2; TEMP = LINK->s[LINK->d1-1]; TEMP1 = LINK->s[LINK->d2-1]; LINK->sigmad = sqrt(TEMP * TEMP / LINK->n[LINK->d1-1] + TEMP1 * TEMP1 / LINK->n[LINK->d2-1]); LINK->t = fabs(LINK->m[LINK->d1-1] - LINK->m[LINK->d2-1]) / LINK->sigmad; LINK->p = 2.0 * (1.0 - ttestprobability(LINK->t, (double)LINK->df)); /* two tailed test!! */ } Local Void computetheenergy(afile, LINK) _TEXT *afile; struct LOC_themain *LINK; { /* compute the energy for d1 minus d2 */ double Z; /* Z score for energy deviation from zero */ double p; /* probability that the energy deviation is zero */ double tol; /* tolerance of p result */ /* compute the energy */ showingerrors = true; memcpy(pmstring, pmstring1, 5L); LINK->r1.estimate = LINK->m[LINK->d1-1]; LINK->r1.error = LINK->s[LINK->d1-1]; LINK->r2.estimate = LINK->m[LINK->d2-1]; LINK->r2.error = LINK->s[LINK->d2-1]; dividerange(LINK->r1, LINK->r2, &LINK->r1overr2); logrange(LINK->r1overr2, baseofnaturallogs, &LINK->logr1overr2); /* energyfactor converts to bits and then to kcal/mole: */ LINK->energyfactor.estimate = 1.0 / (log(2.0) * 2.3424); LINK->energyfactor.error = 0.0; timesrange(LINK->logr1overr2, LINK->energyfactor, &LINK->deltadeltaG); fprintf(afile->f, "deltadeltaG energy, deltaG[%ld] - deltaG[%ld]: $", LINK->d1, LINK->d2); /* write(afile,'deltadeltaG energy, d1 - d2: $');*/ /* writerange(afile,deltadeltaG, round(wid/2), dec, false);*/ writerange(afile, LINK->deltadeltaG, 1L, 2L, false); fprintf(afile->f, "$ kcal/mole\n"); /* write(afile,' r1/r2:');*/ /* writerange(afile,r1overr2, wid, dec, false);*/ /* writeln(afile);*/ /*zzz*/ Z = LINK->deltadeltaG.estimate / LINK->deltadeltaG.error; simpson(Z, &p, &tol); fprintf(afile->f, "Z score for energy deviation from zero: %1.2f\n", Z); fprintf(afile->f, "probability that this deltadeltaG is zero: %1.2f\n", p); /*recompute*/ /* logrange (r1overr2, 2, logr1overr2); (* energyfactor converts to bits and then to kcal/mole: *) energyfactor.estimate := 1.0 /(2.3424); energyfactor.error := 0.0; timesrange(logr1overr2, energyfactor, deltadeltaG); write(afile,'deltadeltaG energy, d1 - d2: $'); writerange(afile,deltadeltaG, round(wid/2), 2, false); writeln(afile,'$ kcal/mole'); */ /*zzz*/ } #undef baseofnaturallogs /* end module testttestprobability */ /******************************************************************************/ /* begin module multtest.themain */ Static Void themain(distributions, multtestp, list, xyin, key) _TEXT *distributions, *multtestp, *list, *xyin, *key; { /* the main procedure of the program */ struct LOC_themain V; Char controlchar; /* character to control the program */ long d = 0; /* number of distributions */ boolean dolog = false; /* do log of means */ double fudgefactor = 1.0; /* fudge factor on standard deviations */ long keyrange = 10; /* range of the key */ _TEXT TEMP; long FORLIM; printf("multtest %4.2f\n", version); if (*list->name != '\0') { if (list->f != NULL) list->f = freopen(list->name, "w", list->f); else list->f = fopen(list->name, "w"); } else { if (list->f != NULL) rewind(list->f); else list->f = tmpfile(); } if (list->f == NULL) _EscIO2(FileNotFound, list->name); SETUPBUF(list->f, Char); fprintf(list->f, "multtest %4.2f\n", version); if (*distributions->name != '\0') { if (distributions->f != NULL) distributions->f = freopen(distributions->name, "r", distributions->f); else distributions->f = fopen(distributions->name, "r"); } else { /* read in the distributions */ rewind(distributions->f); } if (distributions->f == NULL) _EscIO2(FileNotFound, distributions->name); RESETBUF(distributions->f, Char); while (!BUFEOF(distributions->f)) { if (P_peek(distributions->f) == ' ' || P_peek(distributions->f) == '9' || P_peek(distributions->f) == '8' || P_peek(distributions->f) == '7' || P_peek(distributions->f) == '6' || P_peek(distributions->f) == '5' || P_peek(distributions->f) == '4' || P_peek(distributions->f) == '3' || P_peek(distributions->f) == '2' || P_peek(distributions->f) == '1' || P_peek(distributions->f) == '0') { d++; if (d > maxdistrib) { printf("cannot have more than %ld distributions\n", (long)maxdistrib); halt(); } fscanf(distributions->f, "%ld%lg%lg%*[^\n]", &V.n[d-1], &V.m[d-1], &V.s[d-1]); getc(distributions->f); V.s[d-1] = fudgefactor * V.s[d-1]; /* fudge factor for uncertainties */ continue; } if (P_peek(distributions->f) == '*') { fscanf(distributions->f, "%*[^\n]"); getc(distributions->f); continue; } if (P_peek(distributions->f) != 'f') { printf("distributions can only contain numbers or f(udge)\n"); halt(); } getc(distributions->f); fscanf(distributions->f, "%lg%*[^\n]", &fudgefactor); getc(distributions->f); printf("Fudge factor is set to: %*.*f\n", wid, dec, fudgefactor); /* set fudge factor */ } V.colortype = 'n'; V.intervals = 20; if (*multtestp->name != '\0') { if (multtestp->f != NULL) multtestp->f = freopen(multtestp->name, "r", multtestp->f); else multtestp->f = fopen(multtestp->name, "r"); } else rewind(multtestp->f); if (multtestp->f == NULL) _EscIO2(FileNotFound, multtestp->name); RESETBUF(multtestp->f, Char); while (!BUFEOF(multtestp->f)) { if (P_peek(multtestp->f) == ' ' || P_peek(multtestp->f) == '9' || P_peek(multtestp->f) == '8' || P_peek(multtestp->f) == '7' || P_peek(multtestp->f) == '6' || P_peek(multtestp->f) == '5' || P_peek(multtestp->f) == '4' || P_peek(multtestp->f) == '3' || P_peek(multtestp->f) == '2' || P_peek(multtestp->f) == '1' || P_peek(multtestp->f) == '0') { fscanf(multtestp->f, "%ld%ld%*[^\n]", &V.d1, &V.d2); getc(multtestp->f); if (V.d1 > d || V.d2 > d || V.d1 < 1 || V.d2 < 1) { printf( "d1 (=%ld) or d2 (=%ld) is out of the range of distributions, 1 to %ld\n", V.d1, V.d2, d); halt(); } compute(&V); /* writeln(output, 'compare ',d1:1,' to ', d2:1); */ printf("d1 = %2ld: %*ld %*.*f %*.*f\n", V.d1, wid, V.n[V.d1-1], wid, dec, V.m[V.d1-1], wid, dec, V.s[V.d1-1]); printf("d2 = %2ld: %*ld %*.*f %*.*f\n", V.d2, wid, V.n[V.d2-1], wid, dec, V.m[V.d2-1], wid, dec, V.s[V.d2-1]); printf("t = %*.*f\n", wid, dec, V.t); printf("p(%ld=%ld) = %*.*f\n", V.d1, V.d2, wid, dec, V.p); fprintf(list->f, " "); fprintf(list->f, "distribution %ld", V.d1); fprintf(list->f, " | "); fprintf(list->f, "distribution %ld\n", V.d2); fprintf(list->f, "number "); fprintf(list->f, "%*ld", wid, V.n[V.d1-1]); fprintf(list->f, " | "); fprintf(list->f, "%*ld\n", wid, V.n[V.d2-1]); fprintf(list->f, "mean "); fprintf(list->f, "%*.*f", wid, dec, V.m[V.d1-1]); fprintf(list->f, " | "); fprintf(list->f, "%*.*f\n", wid, dec, V.m[V.d2-1]); fprintf(list->f, "standard dev. "); fprintf(list->f, "%*.*f", wid, dec, V.s[V.d1-1]); fprintf(list->f, " | "); fprintf(list->f, "%*.*f\n", wid, dec, V.s[V.d2-1]); fprintf(list->f, "sigma-D = %*.*f\n", wid, dec, V.sigmad); fprintf(list->f, "degrees of freedom = %ld\n", V.df); fprintf(list->f, "t = %*.*f\n", wid, dec, V.t); fprintf(list->f, "p(%ld=%ld) = %*.*f\n", V.d1, V.d2, wid, dec, V.p); fprintf(list->f, "This is the probability that they are the same.\n"); TEMP.f = stdout; *TEMP.name = '\0'; computetheenergy(&TEMP, &V); printf("------------------------------------\n"); computetheenergy(list, &V); fprintf(list->f, "------------------------------------\n\n"); continue; } controlchar = getc(multtestp->f); if (controlchar == '\n') controlchar = ' '; if (controlchar != 'c') { if (controlchar == 'l') dolog = true; else { printf("unknown control character in multtestp: %c\n", controlchar); halt(); } continue; } V.colortype = getc(multtestp->f); if (V.colortype == '\n') V.colortype = ' '; if (V.colortype != 'f' && V.colortype != 's' && V.colortype != 'i' && V.colortype != 'n') { printf("colortype in multtestp must be one of: nisf\n"); halt(); } if (V.colortype == 'i') { fscanf(multtestp->f, "%ld%*[^\n]", &V.intervals); getc(multtestp->f); continue; } if (V.colortype == 'f') { fscanf(multtestp->f, "%lg%*[^\n]", &fudgefactor); getc(multtestp->f); } else { fscanf(multtestp->f, "%*[^\n]"); getc(multtestp->f); } } /* convert to logs if needed */ if (dolog) { printf("USING LOGS\n"); FORLIM = d; for (d = 0; d < FORLIM; d++) { printf("d= %2ld n[d] = %ld m[d] = %*.*f s[d] = %*.*f", d + 1, V.n[d], wid, dec, V.m[d], wid, dec, V.s[d]); V.s[d] /= V.m[d]; /* m[d] := ln(m[d])/ln(2); */ V.m[d] = log(V.m[d]); printf(" m'[d] = %*.*f s'[d] = %*.*f\n", wid, dec, V.m[d], wid, dec, V.s[d]); } } else { /* (* begin module logrange *) procedure logrange(a: range; base: real; var loga: range); (* take log base of a as loga *) begin loga.estimate := ln(a.estimate)/ln(base); loga.error := a.error/ a.estimate ; end; (* end module logrange *) */ printf("NOT USING LOGS\n"); } 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) { /* create color graph */ rewind(xyin->f); } else xyin->f = tmpfile(); } if (xyin->f == NULL) _EscIO2(FileNotFound, xyin->name); SETUPBUF(xyin->f, Char); fprintf(xyin->f, "* multtest %4.2f\n", version); fprintf(xyin->f, "* color graph \n"); V.saturation = 0.5; V.brightness = 1.0; for (V.d1 = 1; V.d1 <= d; V.d1++) { for (V.d2 = 1; V.d2 <= d; V.d2++) { compute(&V); setcolor(&V); showcolor(xyin, (double)V.d1, (double)V.d2, V.symbol, &V); } } /* create color key */ if (*key->name != '\0') { if (key->f != NULL) key->f = freopen(key->name, "w", key->f); else key->f = fopen(key->name, "w"); } else { if (key->f != NULL) rewind(key->f); else key->f = tmpfile(); } if (key->f == NULL) _EscIO2(FileNotFound, key->name); SETUPBUF(key->f, Char); fprintf(key->f, "* multtest %4.2f\n", version); fprintf(key->f, "* color key\n"); V.d2 = 0; /* start at 0.05 so that there is a white space to put the empty box */ /*zzz*/ V.p = 0.0; setcolor(&V); showcolor(key, 0.4, 0.5, 'b', &V); V.symbol = 'c'; for (V.d1 = keyrange * 6; V.d1 <= keysmooth * keyrange; V.d1++) { V.p = (double)V.d1 / (keysmooth * keyrange); setcolor(&V); showcolor(key, (double)V.d1 / keysmooth / 10 - colorshift, 0.5, V.symbol, &V); } /* OLD: (* create color key *) rewrite(key); writeln(key,'* multtest ',version:4:2); writeln(key,'* color key'); d2 := 0; keyrange := 10; for d1 := 0 to keyrange do begin p := d1/keyrange; setcolor; showcolor(key,d1/10-colorshift,0.5,symbol); end; */ /* (* create color key, double range *) d1 := 0; keyrange := 20; for d2 := 0 to keyrange do begin p := d2/keyrange; setcolor; showcolor(key,d1,d2,symbol); end; */ /* hand tests t := 63.7; t := 31.8; t := 12.7; t := 3.8; p := ttestprob(df,t); writeln(list,' t = ',t:wid:dec); writeln(list,' p = ',p:wid:dec); writeln(list,'1-p = ',(1-p):wid:3); */ } #undef color1 #undef color2 #undef color3 #undef color4 #undef color5 #undef colorshift #undef keysmooth #undef dec #undef wid #undef maxdistrib /* end module multtest.themain */ main(argc, argv) int argc; Char *argv[]; { PASCAL_MAIN(argc, argv); if (setjmp(_JL1)) goto _L1; key.f = NULL; strcpy(key.name, "key"); xyin.f = NULL; strcpy(xyin.name, "xyin"); list.f = NULL; strcpy(list.name, "list"); multtestp.f = NULL; strcpy(multtestp.name, "multtestp"); distributions.f = NULL; strcpy(distributions.name, "distributions"); /* testttestprobability; */ themain(&distributions, &multtestp, &list, &xyin, &key); _L1: if (distributions.f != NULL) fclose(distributions.f); if (multtestp.f != NULL) fclose(multtestp.f); if (list.f != NULL) fclose(list.f); if (xyin.f != NULL) fclose(xyin.f); if (key.f != NULL) fclose(key.f); exit(EXIT_SUCCESS); } /* End. */