program gaussint(gout, output); (* gaussint: creates a table of values representing the integral of the gaussian distribution from a z value to infinity. Mark Shaner NCI/FCRDC Bldg 469. Room 144 P.O. Box B Frederick, MD 21702-1201 (301) 846-5581 (-5532 for messages) network address: shaner@ncifcrf.gov National Cancer Institute Laboratory of Mathematical Biology (note that this is an old address) 1992 Maintained by: 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/ *) label 1; (* end of program *) const (* begin module version *) version = 1.70; (* of gaussint.p 2007 Mar 21 2007 Mar 21, 1.70: correct module structure (TDS) 1992 Jun 29, 1.69: functional version 1991 Dec 12, 1.00: origin *) (* end module version *) (* begin module describe.integrate *) (* name gaussint: creates a table of values representing the integral of the gaussian distribution from a z value to infinity. synopsis gaussint(gout : out; output: out) files gout: the output of this program, gout contains a table of values of the integral of the gaussian distribution from varying z values to infinity. description This program finds the integral under the gaussian curve using Simpson's method of numerical integration with an end correction and create a table of many of these values. examples documentation for information about this program, see: Miller, Alan R. PASCAL PROGRAMS FOR SCIENTISTS AND ENGINEERS. SYBEX Inc., 1981. see also author Mark Christopher Shaner bugs there are no known bugs technical notes *) (* end module describe.integrate *) var gout: text; (* output file *) (* 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 intgauss *) procedure intgauss(lower, upper: real; var answer, tol: real); (* Perform a numerical integration of the normal curve 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 *) var deltax, (* the distance between each x value and the subsequent one used in the calculations *) endcor, (* the value of the end correction, it is determined using dy/dx *) endsum, (* the sum of the area under the first and last parabolas *) evensum: real; (* the sum of the area under each of the even numbered parabolas *) i: integer; (* a counter *) oddsum, (* the sum of the area under each of the odd numbered parabolas *) pi: real; (* the value of pi *) pieces: integer; (* the number of parabolas under the curve *) prevsum, (* 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 *) temp, (* a temporary variable used for swapping upper and lower when necessary *) sum: real; (* the value of the area under the curve *) 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. *) x: real; (* the independent variable for calculating the value of functions *) begin sum := 0.0; pi := 4.0*arctan(1.0); val := 1/sqrt (2*pi); tol := 1/maxint; if upper < lower then begin (* switch their values *) temp := lower; lower := upper; upper := temp end; pieces := 2; deltax := (upper-lower)/pieces; x := lower + deltax; oddsum := val*exp(-0.5*x*x); evensum := 0.0; endsum := val*exp(-0.5*lower*lower) + val*exp(-0.5*upper*upper); endcor := -(lower*val)*exp(-0.5*lower*lower) - (-(upper*val)* exp(-0.5*upper*upper)); sum := (endsum + 4.0*oddsum)*deltax/3.0; repeat pieces := pieces*2; prevsum := sum; deltax := (upper - lower)/pieces; evensum := evensum + oddsum; oddsum := 0.0; for i := 1 to (pieces div 2) do begin x := lower + deltax*(2.0*i-1.0); oddsum := oddsum + val*exp(-0.5*x*x) end; sum := (7.0*endsum + 14.0*evensum + 16.0*oddsum + endcor*deltax)* deltax/15.0; until abs(sum-prevsum) <= abs(tol*sum); answer := sum; end; (* end module intgauss *) (* begin module probgauss *) procedure probgauss(zvalue: real; var onetail, tol: real); (* Find the probability associated with the integral from a z value (zvalue) to infinity, (a one tail integral) where onetail is the answer, and tol is the tolerance of the calculations and therefore the possible error in the numerical integration. This procedure requires the intgauss procedure. Written by Mark Shaner, 1992. *) begin intgauss(0,zvalue,onetail,tol); onetail := 0.5 - onetail end; (* end module probgauss *) (* begin module gausstest *) procedure gausstest(var gout: text); var tol, onetail : real; upper: real; i,j: integer; (* counters *) begin rewrite(gout); writeln(gout,'gaussint ',version:4:2); writeln(output,'gaussint ',version:4:2); writeln(gout); writeln(gout,'Values of One-Tail integration of the Normal Distribution'); writeln(gout,'Integral is from z to infinity'); writeln(gout); writeln(gout); write(gout,' z | 0'); for i := 1 to 9 do write(gout,i:7); writeln(gout); for i := 1 to 74 do write(gout,'-'); writeln(gout); i := 300; repeat write(gout,i/100:3:1,'| '); for j := 0 to 9 do begin upper := (i+j)/100; probgauss(upper, onetail, tol); write(gout,onetail:6:4,' '); end; writeln(gout); i := i - 10 until i < 0; end; (* end module gausstest *) begin gausstest(gout); writeln(output,'done'); 1: end.