home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / STAT.ZIP / STAT.DOC < prev    next >
Encoding:
Text File  |  1990-08-08  |  40.5 KB  |  1,347 lines
  1.  
  2.  
  3.  
  4.  
  5.  
  6.                        Libraries of Elementary 
  7.                      Statistical and Mathematical
  8.                             Subprograms
  9.                                 for
  10.                       Turbo Pascal Version 5.X
  11.  
  12.  
  13.                             Version 1.0
  14.  
  15.  
  16.  
  17.  
  18.  
  19.  
  20.  
  21.  
  22.                                by
  23.  
  24.                          Norton Associates
  25.  
  26.  
  27.  
  28.  
  29.  
  30.  
  31.  
  32.  
  33.  
  34.  
  35.  
  36.  
  37.  
  38.  
  39.  
  40.  
  41.  
  42.  
  43.  
  44.                         506 Post Oak Dr.
  45.                         Baytown Tx. 77520
  46.  
  47.                       Compuserve 70441,1337
  48. SUMMARY OF FEATURES
  49.  
  50. Version 1.0 of this software contains 30 elementary statistical and 19 
  51. mathematical subprograms.  The MATH Unit was designed to add the 
  52. FORTRAN-77 intrinsic math functions, except for complex math, that are 
  53. not provided by Turbo Pascal Version 5.X.  However, some additional 
  54. math functions have been added to the MATH Unit as appropriate.
  55.  
  56. STAT Unit
  57.  
  58.            - Uniform and normal random number generators
  59.            - 4 different types mean
  60.            - Summary Statistics ( first four moments, etc ) 
  61.            - Percentiles and quartiles
  62.            - Standard errors of parameters for normal distributions
  63.            - Cumulative frequency normal distribution statistics
  64.            - Integral frequency normal distribution statistics
  65.            - Correlation ( product moment and auto lag)
  66.            - Regression ( polynomial and multiple variable )
  67.            - Smoothing ( binomial, moving average and curvilinear )
  68.            - Sorting ( insert and quick )
  69.            - Other tools ( determinant and remove average )
  70.  
  71. MATH Unit
  72.  
  73.            - Radian to degree conversion
  74.            - Normal trigonometric functions ( tan, secant etc )
  75.            - Inverse trigonometric functions  ( atan2, arcsin etc )
  76.            - Hyperbolic functions ( sinh, cosh, tanh )
  77.            - Logarithmic functions ( power, log to any base etc )
  78.            - Other functions (factorial and root of an equation )
  79.  
  80. POSSIBLE FUTURE ENHANCEMENTS
  81.  
  82.            - Cross lag correlation ( two variable)
  83.            - Increase variety of standard errors
  84.            - Contour plotting program
  85.            - Spline smoothing functions
  86.            - Integration
  87.            - Error functions
  88.            - Interpolation and extrapolation
  89.            - Gamma and Bessel functions
  90.            - Graphics ( plot function )
  91.            - Better use of 8087 coprocessor ( i.e., use 8087 tan 
  92.                 function
  93.            - Arrays larger than 16380 single precision elements.  Have 
  94.              created a single dimension array containing 130000 
  95.              elements not 16380 as currently implemented.
  96.            - Complex Math functions
  97. 1.0 Introduction
  98.  
  99. 1.1 Background
  100.  
  101. The statistical programs contained in the STAT Unit were originally 
  102. developed in FORTRAN over the years.  With the advent of a good Pascal 
  103. environment, the use and need of FORTRAN to statistically crunch data 
  104. have diminished drastically. 
  105.   
  106. The MATH Unit was developed because the language of Pascal does not 
  107. contain all the intrinsic math functions as those contained in the 
  108. standard FORTRAN-77 intrinsic function library [6].  The combination 
  109. of Turbo Pascal's math functions and those contained in this MATH Unit 
  110. provides the same functionality as those contained in the FORTRAN-77 
  111. intrinsic function library, except for complex math.
  112.  
  113. The user of these libraries are assumed to be familiar with 
  114. statistics.  The definitions of terms used in these libraries can be 
  115. found by reviewing the reference contained in Section 5.0.  Additional 
  116. statistical and mathematical tools will be added in the future.  
  117.  
  118. This is a shareware product.  Register your copy of the STAT and MATH 
  119. libraries by reading the README.DOC file contained with this 
  120. disk(ette) and sending in the nominal $35 to register the software.  
  121.  
  122.  
  123. 1.2 Converting the Statistical Unit to Turbo Pascal
  124.  
  125. Several things happened while converting the statistical software from 
  126. FORTRAN to Pascal.  
  127.  
  128. First, the size of the arrays - With the the advent of Personal 
  129. Computers (PC's), many mainframe FORTRAN programs have been converted 
  130. to PC's.  In most cases special care had to be taken when the size of 
  131. the arrays multiplied by the size of the element exceeded 64k (65520 
  132. bytes).  This limitation was due to the 8086 segmented architecture. 
  133. Most FORTRAN vendors implemented a way around the 8086 segmented 
  134. architecture.  Well, Turbo Pascal has no provided such a capability 
  135. [1].  The segmented architecture of the 8086 family basically limits 
  136. the size of an array.  The below table explains the size of the 
  137. element and the number of elements in the 65k 8086 architecture 
  138. boundary.
  139.  
  140.                        Table 1. 
  141.  
  142. Element     Size      Maximum number of       
  143. type       (bytes)    elements in array       STAT types
  144. -------    --------   ------------------  ----------------------
  145. Single        4         16380             single_array_pointer
  146. Longint       4         16380             longint_array_pointer
  147.  
  148. For example, this table states that the maximum number of elements 
  149. (indices of the array) of a single precision array is 16380 elements 
  150. in a 65K segment.  This is satisfactory for most applications.  
  151. However, for some applications a vast increase in the number of 
  152. elements in the array was needed.  To make a long story short, I have 
  153. developed a single precision array which which contains 130000 
  154. elements in the array.  This implies that the memory requirements for 
  155. this array size alone is 520000 bytes.  I have not completed testing 
  156. of this version of the software.
  157.  
  158. Second, the 640k DOS limitation - There are two major issues.  How can 
  159. a user make best use of the 640k DOS limitation and if the user needed 
  160. more memory could he/she use extended/expanded memory for array 
  161. storage?  The latter issue is not addressed in this version of the 
  162. software.  The former issue was challenging enough.
  163.  
  164. With a fixed size of DOS, the number of arrays and the number of 
  165. elements in the arrays is very finite.  Not knowing the number of 
  166. elements required by each user dictated that the STAT unit use the 
  167. maximum number of elements from Table 1.  If the arrays were 
  168. statically created this meant that only 9 arrays of 65k bytes could 
  169. ever be created.  This was totally not satisfactory.  The solution was 
  170. to dynamically create the arrays on the heap.  The user could then 
  171. have the best of both worlds.  Basically, heap is the part of memory 
  172. between the top of the application in memory and the 640k DOS 
  173. limitation.  By dynamically creating the arrays, the number of arrays 
  174. in the application is now very large and the size of each array may be 
  175. up to the full 65k limit.  More information on this subject is 
  176. contained in the next section.
  177.  
  178. Finally, the precision - To make the STAT unit more precise, all 
  179. normally double precision numbers are in extended precision.  By 
  180. adding two more bytes per variable, the number of significant digits 
  181. increased by four (15 to 19 digits).
  182.  
  183.  
  184. 1.3 User's Guide to the STAT and MATH Units
  185.  
  186.  
  187. 1.3.1 Format for the Libraries in Sections 3.0 and 4.0
  188.  
  189. A standard format for presenting the subprograms contained in the 
  190. libraries have been developed for Sections 3.0 and 4.0.  An outline of 
  191. the standard format is presented directly below.
  192.  
  193. Name of subprogram here
  194.  
  195. Purpose:   A short purpose statement goes here
  196.  
  197. Interface: The Turbo Pascal Unit Interface
  198.  
  199. Errors:    Any errors are reported here
  200.  
  201. Remarks:   Any special considerations go here
  202.  
  203. Example:   Sample code go here
  204.  
  205.  
  206. 1.3.2  STAT Unit
  207.  
  208. To prevent any unforeseen problems, a series of STAT Unit procedures 
  209. have been developed to dynamically create and delete arrays.  These 
  210. procedures are listed below.
  211.  
  212.  
  213.                     Table 2.
  214.  
  215.                           Procedures for 
  216. Element      Dynamic Creation and Deletion Procedures
  217. Type               Create              Delete
  218. ---------    -------------------  -------------------
  219. Single       create_single_array  delete_single_array
  220. Longint     create_longint_array  delete_longint_array
  221.  
  222.  
  223. The creation and deletion of arrays is shown in the demo called 
  224. DYNAMIC.PAS.  The user should execute DYNAMIC.EXE.  Basically, the 
  225. user must initiate and create the array before any information is put 
  226. into it.  To initiate an array the user must know what type of an 
  227. array you want.  For example, if the user wants a single dimensioned 
  228. array, he must initiate the array in the VAR section of the Pascal 
  229. code with one of the STAT types mentioned in Table 1.  Once Initiated 
  230. the user may then create the array.  This is done by selecting 
  231. corresponding creation procedures from Table 2.  Once created, the 
  232. user may then use the array as any other array.  Because the arrays 
  233. are created on the heap, the identification of the array is not 
  234. "lead[1]" but "lead^[1]".  The caret (^) indicates to Pascal that the 
  235. array was created on the heap.  The following example demonstrates the 
  236. creation process. 
  237.  
  238.  
  239. program demo;
  240.  
  241. const
  242.      max = 15000;                   { 15000 elements in the array }
  243. var
  244.      lead : single_array_pointer;   { Step 1: initiate a single type }
  245. begin
  246.      create_single_array(max,lead); { Step 2: create the lead array of 
  247.                                          max elements }
  248.      for j := 1 to max do           { Step 3: use the array normally }
  249.            lead^[j] := j;
  250. end. 
  251.  
  252. The Pascal programs MULREG.PAS and TIMER.PAS are additional demo 
  253. programs to see the STAT Unit in action.  The Pascal constants and 
  254. types for the MATH and STAT Units are defined in Section 2.0.  All 
  255. subprograms of the STAT Unit are discussed using the format described 
  256. above are in Section 3.0.  
  257.  
  258.  
  259. 1.3.3 MATH Unit
  260.  
  261. The MATH Unit functions are designed to perform just like the math 
  262. intrinsic function contained in FORTRAN-77.  The Pascal program 
  263. MATHTEST.PAS is provided as a demo to test most MATH unit subprograms.  
  264. The MATH Unit constants are defined in Section 2.0 and a description 
  265. of these functions are contained in Section 4.0 
  266. 2.0 Predefined Constants and Types
  267.  
  268.  
  269. 2.1 Predefined Library Constants
  270.  
  271. The following constants are defined in the STAT and MATH units.  The 
  272. user may use these constants as needed in programs.  Make sure that 
  273. you do not change these constant values.  Be careful, Turbo Pascal 
  274. does allow the programmer to change the values defined in the constant 
  275. section.
  276.  
  277. STAT Unit
  278.  
  279.            error      = -99.9;                    { for most errors}
  280.            maxsize    = 65520;                    { max segment size }
  281.            maxsingle  = maxsize div (sizeof(single));  { 16380 }
  282.            maxlongint = maxsize div (sizeof(longint)); { 16380 }
  283.            v1         = 1;                       { constants }
  284.            v2         = 10000;                   { for uniform random}
  285.            v3         = 3000;                    {  number generators}
  286.            maxorder   = 10;                 { max order of polynomial}
  287.            maxmatrix  = maxorder * 2 - 1; (19)
  288.            zero       = 0.0;
  289.            one        = 1.0;
  290.            two        = 2.0;
  291.            uno        = 1;
  292.            dos        = 2;
  293.            c2         = 2.0;
  294.            c3         = -3.0;
  295.            c4         = 4.0;
  296.            c5         = 5.0;
  297.            c6         = 6.0;
  298.            c11        = 11.0;
  299.            c12        = 12.0;
  300.            i_4        = 1.0/c4; { 0.25 };
  301.            i_16       = 1.0/16.0;
  302.            i_35       = 1.0/35.0;
  303.  
  304. MATH Unit
  305.  
  306.            pi         = 3.14159265359;
  307.            pi_2       = pi / 2.0; 
  308.            pi2        = pi * 2.0;
  309.            rad        = pi / 180.0
  310.            i_rad      = 180.0 / pi;
  311.            one        = 1.0;
  312.            zero       = 0.0;
  313.            infinity   = 1.0E+09;
  314.  
  315.  
  316. 2.2 Predefined Types
  317.  
  318. STAT Unit
  319.            single_array_dummy    = array[1..maxsingle] of single;
  320. *          single_array_pointer  = ^single_array_dummy;
  321.  
  322.            longint_array_dummy   = array[1..maxlongint] of longint;
  323. *          longint_array_pointer = ^longint_array_dummy;     
  324.  
  325. *          quartype       = array[1..5] of single;
  326. *          arry_type      = array[1..maxorder,1..maxorder] of extended
  327.  
  328. ------------
  329. * Only use these in VAR section of Turbo Pascal program
  330.  
  331.  
  332. MATH Unit 
  333.  
  334.            None
  335. 3.0 Statistical Unit
  336.  
  337. 3.1 Random Number Generators
  338.  
  339.  
  340. 3.1.1 Uniform1
  341.  
  342. Purpose:   Calculate a uniform distribution floating point number 
  343.            between 0.0 and 1.0 [12].
  344.  
  345. Interface: function uniform1: single;
  346.  
  347. Errors:    None
  348.  
  349. Remarks:   v1,v2 and v3 are initial seeds that are defined as 
  350.            constants and may be altered for different random number 
  351.            sequences.  The domain of these seeds are between 1 and 
  352.            32767. 
  353.  
  354. Example:   begin
  355.                  x := uniform1;
  356.            end;
  357.  
  358.  
  359. 3.1.2 Rndnormal1
  360.  
  361. Purpose:   Generate a floating point number whose sample population is 
  362.            normal and between +- 6 standard deviations
  363.  
  364. Interface: function rndnormal1(mean , standev : extended) : single;
  365.            
  366.            Where: mean    = mean of sample population and
  367.                   standev = standard deviation of sample population
  368.  
  369. Errors:    None
  370.  
  371. Remarks:   Rndnormal1 is based on the polar coordinate normal random 
  372.            number generator [11].
  373.  
  374. Example:   begin
  375.                 mean = 100.0;
  376.                 sd   =  5.0;
  377.                 x := rndnormal1(mean,sd);
  378.            end;
  379.  
  380.  
  381. 3.1.3 Rndnormal2
  382.  
  383. Purpose:   Provide an approximate normal distribution random floating 
  384.            point number between +- 6 standard deviations
  385.  
  386. Interface: function rndnormal2(mean , standev : extended) : single;
  387.  
  388.            Where: mean    = mean of expected sample population
  389.                   standev = standard deviation of expected sample 
  390.                              population
  391.  
  392. Errors:    None
  393.  
  394. Remarks:   This normal random number is an approximation to the normal 
  395.            random discussed previously [5].
  396.  
  397. Example:   begin
  398.                 mean = 100.0;
  399.                 sd   = 5.0;
  400.                 x := rndnormal2(mean,sd);
  401.            end;
  402. 3.2 Sorting
  403.  
  404. 3.2.1 Insert Sort
  405.  
  406. Purpose:   Classical data structure for sorting floating point 
  407.            numbers.
  408.  
  409. Interface: procedure insert(n: word; var a: single_array_pointer);
  410.  
  411.            Where: n = number of values to be sorted
  412.                   a = array of values to be sorted
  413.            
  414. Errors:    insert> no sample data - Number of sample data points less 
  415.            than two.
  416.  
  417. Remarks:   For arrays up to approximately 15 floating point numbers 
  418.            the insert sort is more efficient than quicksort.  The code 
  419.            and data requirements of insert are smaller than quicksort.  
  420.            The original order of the input array is destroyed.
  421.  
  422. Example:   begin
  423.                 insert(n,a);
  424.            end;
  425.  
  426.  
  427. 3.2.2 Quick sort 
  428.  
  429. Purpose:   Provided a very fast way for sorting floating point 
  430.            numbers.
  431.  
  432. Interface: procedure qsort(n: word; var a: single_array_pointer);
  433.  
  434.            Where: n = number of values to be sorted
  435.                   a = array of values to be sorted
  436.            
  437. Errors:    1) quick> no sample data - Number of sample data point less 
  438.            than 2.
  439.  
  440.            2) quick> stack space exceeded - number of stacks used in 
  441.            quicksort has been exceeded. 
  442.  
  443. Remarks:   Quick sort provides a very fast way of sorting floating 
  444.            point numbers.  The original order of the input array is 
  445.            destroyed.
  446.  
  447. Example:   begin
  448.                 quick(n,a);
  449.            end;
  450. 3.3  Means
  451.  
  452. 3.3.1 Means
  453.  
  454. Purpose:   Calculate the arithmetic, geometric, harmonic and root mean 
  455.            square means for the sample data [8].
  456.  
  457. Interface: procedure means(n: word; a : single_array_pointer;  
  458.                           var xmean, gmean, hmean, rmsmean : single);
  459.  
  460.            Where: xmean  = arithmetic mean
  461.                   gmean  = geometric mean
  462.                   hmean  = harmonic mean
  463.                  rmsmean = root mean square mean
  464.  
  465. Errors:    1) means> no sample data: n < 2 - Less than two data points 
  466.            provided.
  467.  
  468.            2) means> Gmean can not be calculated - There is a zero or 
  469.            negative number in the sample data.
  470.  
  471.            3) means> Hmean can not be calculated - There is a zero or 
  472.            negative number in the sample data
  473.  
  474. Remarks:   If there are no variations in the data then all means have 
  475.            the same value.  As the variation in the data increases, 
  476.            the difference between the means increase.  Each type of 
  477.            mean have their own advantages and disadvantages.
  478.  
  479. Example:   begin
  480.                 means(n,a,xmean,gmean,hmean,rmsmean);
  481.            end;
  482.  
  483. 3.4 Summary statistics
  484.  
  485.  
  486. 3.4.1 Elementary statistics
  487.  
  488. Purpose:   Calculate elementary statistics of data
  489.  
  490. Interface: elem_stat(n : word; a : single_pointer_array;
  491.                 var small, large, mean, sd : single);
  492.  
  493.            Where: small = smallest value 
  494.                   large = largest value
  495.                   mean  = arithmetic mean
  496.                   sd    = adjusted standard deviation for sample 
  497.                           population (i.e., n - 1).
  498.  
  499. Errors:    1) elem_stat> no sample data - 0 or negative number of data 
  500.            points.
  501.  
  502.            2) elem_stat> number of sample = 1: no standard deviation - 
  503.            standard deviation can not be calculated with one data 
  504.            point. 
  505.  
  506. Remarks:   Elem_stat is a very simple way of calculating the simple 
  507.            statistics of the sample population.  Other statistical 
  508.            parameters may be calculated from this data (i.e., range = 
  509.            large - small + 1).   
  510.  
  511. Example:   begin
  512.                 elem_stat(n,a,small,large,mean,sd);
  513.            end;
  514.  
  515.  
  516. 3.4.2 Weighted arithmetic mean
  517.  
  518. Purpose:   Calculate a frequency or weighted summary statistics
  519.  
  520. Interface: wxmean(n : word; a : single_array_pointer; freq : 
  521.                 longint_array_pointer; var small, large, mean, sd : 
  522.                 single);
  523.  
  524.            Where: freq = interval count of data for each corresponding 
  525.                           value in the a array
  526.  
  527. Errors:    1) wxmean> no sample data point - number of data points 
  528.            less than one.
  529.  
  530.            2) wxmean> no standard deviation calculated - Standard 
  531.            deviation can not be calculated from one data point. 
  532.            
  533. Remarks:   Do not forget that freq contains the frequencies for each 
  534.            sample value in the a array, respectively.
  535.  
  536. Example:   begin
  537.                 wxmean(n,a,freq,small,large,mean,sd);
  538.            end;
  539.  
  540.  
  541. 3.4.3 Moments
  542.  
  543. Purpose:   Calculate the first 4 moments for a sample population
  544.  
  545. Interface: procedure moments( n:word; a : single_array_pointer;
  546.                           var mean,sd,skew,beta2: single);
  547.            Where: skew = skewness
  548.                  beta2 = beta2.
  549.  
  550. Errors:    1) moments> no sample data - number of sample data point 
  551.            less than 2.
  552.  
  553.            2) moments> no data for beta2 and skewness - the standard 
  554.            deviation was zero. 
  555.  
  556. Remarks:   This procedure provides a simple way of calculating the 
  557.            first through the fourth moments.  Qurtosis = beta2 - 3;
  558.  
  559. Example:   var a : single_array_pointer;
  560.            begin
  561.                 create_single_array(n,a);
  562.                 ...put data into array a...;
  563.                 moments(n,a,mean,sd,skew,beta2);
  564.            end;
  565. 3.5 Percentiles
  566.  
  567. 3.5.1 Quartiles
  568.  
  569. Purpose:   Calculate the first through third quartile and smallest and 
  570.            largest value of the sample population
  571.  
  572. Interface: procedure quartile( n: word; a : single_array_pointer;
  573.                           var quart : quartype );
  574.            Where quart is an array of 5 elements containing the 3 
  575.            quartiles and the smallest and largest.
  576.  
  577. Errors:    quartile> no sample data - number of sample data point less 
  578.            than 2.   
  579.  
  580. Remarks:   The smallest        = quart[1]      0.0 percentile
  581.            first quartile      = quart[2]     25.0 percentile
  582.            second quartile     = quart[3]     50.0 percentile
  583.            third quartile      = quart[4]     75.0 percentile
  584.            largest             = quart[5]    100.0 percentile
  585.  
  586. Example:   var  quart: quartype:
  587.            begin
  588.                 quartile(n,a,quart);
  589.            end;
  590.  
  591.  
  592. 3.5.2 Percentile
  593.  
  594. Purpose:   Calculate a percentile of the sample population
  595.  
  596. Interface: function percentile( n : word; a : single_array_pointer;
  597.                           percent : single ): single;
  598.  
  599. Errors:    percentile> no sample data - Number of sample data points 
  600.            less than 2.
  601.  
  602. Remarks:   Make sure the variable percent is a percentile (e.g., %)
  603.  
  604. Example:   begin
  605.                 x := percentile(n,a,23.0);
  606.            end;
  607.  
  608. 3.6 Standard error
  609.  
  610. Purpose:   Calculate the standard error [8]. 
  611.  
  612. Interface: function standard_error( n: word, sd : single ;
  613.                                     ntype : word) : single;
  614.  
  615.            Where   ntype  standard error of
  616.                    -----  -----------------------
  617.                      1    mean
  618.                      2    sum 
  619.                      3    median
  620.                      4    standard deviation
  621.                      5    mean deviation
  622.                      6    coefficient of variation
  623.                      7    coefficient of correlation
  624.  
  625. Error:     1) standard_error> type out of range - ntype must be 1 - 7 
  626.            inclusive
  627.  
  628.            2) standard_error> no sample data - number of samples less 
  629.            than 2
  630.  
  631. Remarks:   These standard errors are based on a normal distribution.
  632.  
  633. Example:   begin
  634.                 x := standard_error(n,sd,1);   { for the mean }
  635.            end;
  636.  
  637.  
  638. 3.7 Cumulative Distribution function: Normal distribution
  639.  
  640.  
  641. 3.7.1 From probability to standard deviation
  642.  
  643. Purpose:   Calculate the correct standard deviation from a given 
  644.            probability [14]
  645.  
  646. Interface: function cdf_prob_to_sd(prob: single): single;
  647.  
  648. Errors:    cdf_prob_to_sd> input probability out of range - input 
  649.            probability must between 0 and 1. 
  650.  
  651. Remarks:   The algorithm used is accurate to 4 decimal places
  652.  
  653. Example:   begin
  654.                 x := cdf_prob_to_sd(0.50);   { x = 0.0 }
  655.            end;
  656.  
  657.  
  658. 3.7.2 From standard deviation to probability
  659.  
  660. Purpose:   Calculate the correct probability from a given standard 
  661.            deviation [14]
  662.  
  663. Interface: function cdf_sd_to_prob(sd: single): single;
  664.  
  665. Errors:    None 
  666.  
  667. Remarks:   The algorithm used is accurate to 4 decimal places
  668.  
  669. Example:   begin
  670.                 x := cdf_sd_to_prob(0.0);   { x = 0.50 }
  671.            end;
  672.  
  673. 3.8 Integral Distribution function: Normal distribution
  674.  
  675. 3.8.1 From probability to standard deviation
  676.  
  677. Purpose:   Calculate the correct standard deviation from a given 
  678.            probability [14]
  679.  
  680. Interface: function int_prob_to_sd(prob: single): single;
  681.  
  682. Errors:    int_prob_to_sd> input probability out of range - 
  683.            probability must be between 0 and 1.
  684.  
  685. Remarks:   The algorithm used is accurate to 4 decimal places
  686.  
  687. Example:   begin
  688.                 x := int_prob_to_sd(0.66);   { x = 1.0 }
  689.            end;
  690.  
  691.  
  692. 3.8.2 From standard deviation to probability
  693.  
  694. Purpose:   Calculate the correct probability from a given standard 
  695.            deviation [14]
  696.  
  697. Interface: function int_sd_to_prob(sd: single): single;
  698.  
  699. Errors:    None
  700.  
  701. Remarks:   The algorithm used is accurate to 4 decimal places
  702.  
  703. Example:   begin
  704.                 x := int_sd_to_prob(-1.0);   { x = 0.37 }
  705.            end;
  706.  
  707. 3.9 Correlation
  708.  
  709. 3.9.1 Linear correlation coefficient
  710.  
  711. Purpose:   Calculate the linear correlation coefficients via the 
  712.            product moment technique
  713.  
  714. Interface: procedure corr_coef(n: word; x, y : single_array_pointer; 
  715.                           var r : single);
  716.  
  717.            Where: x and y are the independent and dependent arrays, 
  718.            respectively.
  719.            r = correlation coefficient
  720.  
  721. Errors:    corr_coef> no sample data - number of data point less than 
  722.            2.
  723.  
  724. Remarks:   The correlation coefficient is adjusted for sample size 
  725.            (e.g., n - 1)
  726.  
  727. Example:   begin
  728.                 ... get x and y data...;
  729.                 corr_coef(n,x,y,r);
  730.            end;
  731.  
  732.  
  733. 3.9.2 Linear Auto lag correlation coefficient
  734.  
  735. Purpose:   Calculate the linear correlation coefficient via the 
  736.            product moment technique for each step of the lag. 
  737.  
  738. Interface: procedure auto_corr(n: word; x : single_array_pointer; 
  739.                      lag: word; var auto : single_array_pointer);
  740.  
  741.            Where: lag = number of lags to calculate the auto 
  742.                           correlation coefficients 
  743.                  auto = an array of lag number of elements
  744.  
  745. Errors:    1) auto_corr> no sample data - number of sample data points 
  746.            less than 2.
  747.  
  748.            2) auto_corr> lag can not exceed number in sample data - no 
  749.            resultant data points can be used.
  750.  
  751. Remarks:   The correlation coefficient is adjusted for sample size 
  752.            (e.g., n - 1)
  753.  
  754. Example:   var x,auto : single_array_pointer
  755.            begin
  756.                 ... get x and y data...;
  757.                 auto_corr(n,x,20,auto);
  758.            end;
  759. 3.10 Regression : least squares
  760.  
  761.  
  762. 3.10.1 Linear fit
  763.  
  764. Purpose:   Calculate the linear regression coefficients and 
  765.            correlation coefficient for pairs of x and y values [9].
  766.  
  767. Interface: procedure linfit(x,y,sigmay : single_array_pointer ;
  768.                           n,mode: word; var a,b,r : single);
  769.  
  770.            Where: sigmay = standard deviation of measurement for each 
  771.                                y variable
  772.                  if mode = 0 then sigmay array not used
  773.                  if mode = 1 then sigmay array is used
  774.                  a = intercept
  775.                  b = slope
  776.                  r = linear correlation coefficient
  777.                  y estimated =  a + bx
  778.  
  779. Errors:    linfit> number of data points < 2 - number of data points 
  780.            used in linfit is less than 2
  781.  
  782. Remarks:   Sigmay is not used in linfit if mode = 0
  783.  
  784. Example:   var x,y,sigmay : single_array_pointer;
  785.            begin
  786.                 ... create x and y arrays...;
  787.                 ... put data into x and y arrays...;
  788.                 linfit(x,y,sigmay,n,0,a,b,r);
  789.            end;
  790.  
  791.  
  792. 3.10.2 Polynomial fit
  793.  
  794. Purpose:   Calculate the polynomial regression coefficients and 
  795.            multiple regression correlation coefficient for pairs of x 
  796.            and y values [9].
  797.  
  798. Interface: procedure polfit(n: word ;x, y, sigmay : 
  799.            single_array_pointer ; nterms, mode, : word; var a: 
  800.            single_array_pointer; r,se : single);
  801.  
  802.            Where: sigmay is the same as in linfit
  803.                 nterms = order + 1. For parabola nterms = 3
  804.                 a = array of regression coefficients
  805.                 r = multiple regression correlation coefficient
  806.                 se = standard error of regression
  807.                 y estimated = a[1] + a[2]*x etc
  808. Errors:    1) polfit> number of data point < 2 - number of data points 
  809.            used in polfit is less than 2.
  810.  
  811.            2) polfit> maxorder of polfit is 10 - nterms was larger 
  812.            than 10.  
  813.  
  814.            3) polfit> number of terms exceeds the number of data 
  815.            points.  The equations used in polfit required that the 
  816.            number of data points exceed the number of data points.
  817.  
  818.            4) polfit> determinant or correlation coefficient 
  819.            approximately 0 - Polfit became unstable and reported the 
  820.            error.
  821.  
  822. Remarks:   Sigmay is never used in polfit if mode = 0
  823.  
  824. Example:   var x,y,a : single_array_pointer;
  825.            begin
  826.                 ... create x and y arrays...;
  827.                 ... put data into x and y arrays...;
  828.                 polfit(n,x,y,sigmay,3,0,a,r,se);  ( parabola)
  829.            end;
  830.  
  831.  
  832. 3.10.3 Multiple Regression fit
  833.  
  834. Purpose:   Calculate the multiple regression coefficients and multiple 
  835.            regression correlation coefficient for pairs of x, y and z 
  836.            values [7].
  837.  
  838. Interface: procedure mul_reg(n : word; x, z ,y : single_array_pointer 
  839.                 ; var a: single_array_pointer; r,se : single);
  840.  
  841.            Where: z = one of two arrays of independent measurements
  842.                   a = array of regression coefficients
  843.                   r = multiple regression correlation coefficient
  844.                   se = standard error of regression
  845.                   y estimated = a[1] + a[2]*x + a[3]*z
  846.  
  847. Errors:    1) mul_reg> sample data less than 4 - because this uses 3 
  848.            constraints then number of data points required is 4.
  849.  
  850.            2) mul_reg> determinant or correlation coefficient 
  851.            approximately 0 - Mul_reg became unstable and reported the 
  852.            error.
  853.  
  854. Remarks:   Mul_reg is used for two independent and one dependent 
  855.            variable arrays.
  856.  
  857. Example:   var x,z,y,a : single_array_pointer;
  858.            begin
  859.                 ... create x and y arrays...;
  860.                 ... put data into x and y arrays...;
  861.                 mul_reg(n,x,z,y,a,r,se);
  862.            end;
  863.  
  864. 3.11 Data Smoothing 
  865.  
  866. 3.11.1 Binomial 1-2-1
  867.  
  868. Purpose:   Smooth data with a binomial 1-2-1 function [9]
  869.  
  870. Interface: procedure smooth121(n:word; var a: single_array_pointer);
  871.  
  872. Errors:    smooth121> no sample data - number of data points less than 
  873.            3.
  874.  
  875. Remarks:   This procedure smooths out aberrations in the data
  876.  
  877. Example:   begin
  878.                 smooth121(n,a);
  879.            end;
  880.  
  881.  
  882. 3.11.2 Binomial 1-4-6-4-1
  883.  
  884. Purpose:   Smooth data with a binomial 1-4-6-4-1 function [9]
  885.  
  886.  
  887. Interface: procedure smooth14641(n:word; var a: single_array_pointer);
  888.  
  889. Errors:    smooth14641> no sample data - number of data points less 
  890.            than 5.
  891.  
  892. Remarks:   This procedure has approximately twice the smoothing power 
  893.            of smooth121.
  894.  
  895. Example:   begin
  896.                 smooth14641(n,a);
  897.            end;
  898.  
  899.  
  900. 3.11.3 Curvilinear
  901.  
  902. Purpose:   Smooth data which is markedly concave or convex [13]
  903.  
  904. Interface: procedure smoothcurve(n:word; var a: single_array_pointer);
  905.  
  906. Errors:    smoothcurve> no sample data - number of data points less 
  907.            than 6.
  908.  
  909. Remarks:   This procedure smooths data points that are from a very 
  910.            concave or convex function.
  911.  
  912. Example:   begin
  913.                 smoothcurve(n,a);
  914.            end;
  915.  
  916.  
  917. 3.11.4 Moving Average
  918.  
  919. Purpose:   Smooth data using moving average approach
  920.  
  921.  
  922. Interface: function mov_avg(n:word; var a: single_array_pointer; ma, 
  923.                      counter: word);
  924.            Where: ma = number of elements in moving average.
  925.              counter = current pointer in array
  926.  
  927. Errors:    movavg> number of data points less than moving average - 
  928.            The starting point of the calculation is less than the 
  929.            first data point.
  930.  
  931. Remarks:   None
  932.  
  933. Example:   begin
  934.                 x := mov_avg(n,a,ma,j);
  935.            end;
  936. 3.12 Other tools
  937.  
  938.  
  939. 3.12.1 Determinant
  940.  
  941. Purpose:   Calculate the determinant [3]
  942.  
  943.  
  944. Interface: function determ(arry : arry_type ; norder : integer) : 
  945.                      extended
  946.            Where: arry = array of values to calculate the determinate.
  947.                 norder = order of the arry
  948.  
  949. Errors:    determinate> Determinate contains a zero - Can not solve 
  950.            equation with a zero.
  951.  
  952. Remarks:   None
  953.  
  954. Example:   const     norder = 3;
  955.            var       arry : arry_type;
  956.            begin
  957.                 stuff data into arry   
  958.                 x := determ(arry,norder);
  959.            end;
  960.  
  961. 3.12.2 Remove Average
  962.  
  963. Purpose:   Remove the average or any constant from the values in the 
  964.            data array
  965.  
  966. Interface: function remove_avg(n:word; var a: single_array_pointer;  
  967.                      avg: single);
  968.            Where: avg = value to be deleted from the array.
  969.  
  970. Errors:    Remove_avg> number of data points less than 2 - There to 
  971.            few data points.
  972.  
  973. Remarks:   None
  974.  
  975. Example:   begin
  976.                 x := remove_avg(n,a,avg);
  977.            end;
  978. 4.0 Mathematical
  979.  
  980. 4.1 Radian Conversion
  981.  
  982. 4.1.1 Degrees to Radians
  983.  
  984. Purpose:   To convert from degrees to radian
  985.  
  986. Interface: function deg_rad( x : single) : single);
  987.  
  988. Errors:    None
  989.  
  990. Remarks:   None
  991.  
  992. Example:   begin
  993.                 y := deg_rad(180.0);   (y = 3.14 radians)
  994.            end;
  995.  
  996.  
  997. 4.1.2 Radians to Degrees
  998.  
  999. Purpose:   To convert from radians to degrees
  1000.  
  1001. Interface: function rad_deg( x : single) : single);
  1002.  
  1003. Errors:    None
  1004.  
  1005. Remarks:   None
  1006.  
  1007. Example:   begin
  1008.                 y := rad_deg(pi);   (y = 180.0 degrees)
  1009.            end;
  1010.  
  1011. 4.2 Trigonometric Functions
  1012.  
  1013. 4.2.1 Tan
  1014.  
  1015. Purpose:   Turbo Pascal does not have a Tan
  1016.  
  1017. Interface: function tan( x : single) : single;
  1018.  Errors:   Cosine term can not be zero
  1019.  
  1020. Remarks:   None
  1021.  
  1022. Example:   begin
  1023.                  
  1024.                 x = tan(deg_rad(45.0));  ( x = 1.0)
  1025.            end;
  1026.  
  1027.  
  1028. 4.2.2 Secant
  1029.  
  1030. Purpose:   Turbo Pascal does not have a Secant
  1031.  
  1032. Interface: function secant( x : single) : single;
  1033.  
  1034. Errors:    Cosine term can not be zero
  1035.  
  1036. Remarks:   None
  1037.  
  1038. Example:   begin
  1039.                 x = secant(deg_rad(45.0));  ( x = 1.414)
  1040.            end;
  1041.  
  1042.  
  1043. 4.2.3 Cosecant
  1044.  
  1045. Purpose:   Turbo Pascal does not have a Cosecant
  1046.  
  1047. Interface: function cosecant( x : single) : single;
  1048.  
  1049. Errors:    Sine term can not be zero
  1050.  
  1051. Remarks:   None
  1052.  
  1053. Example:   begin
  1054.                 x = cosecant(deg_rad(45.0));  ( x = 1.414)
  1055.            end;
  1056.  
  1057.  
  1058. 4.2.4 Cotan
  1059.  
  1060. Purpose:   Turbo Pascal does not have a Cotan
  1061.  
  1062. Interface: function cotan( x : single) : single;
  1063.  
  1064. Errors:    Tan term can not be zero
  1065.  
  1066. Remarks:   None
  1067.  
  1068. Example:   begin
  1069.                 x = cotan(deg_rad(45.0));  ( x = 1.0)
  1070.            end;
  1071. 4.3 Inverse Trigonometric Functions
  1072.  
  1073. 4.3.1 Arc sin
  1074.  
  1075. Purpose:   Turbo Pascal does not have a arcsin
  1076.  
  1077. Interface: function arcsin( x : single) : single;
  1078.  
  1079. Errors:    Input value x must be within plus and minus 1.0
  1080.  
  1081. Remarks:   Uses arctan identity to solve for arcsin
  1082.  
  1083. Example:   begin
  1084.                 x = rad_deg(arcsin(sin(deg_rad(45.0))));  ( x = 45.0)
  1085.            end;
  1086.  
  1087.  
  1088. 4.3.2 Arc cosine
  1089.  
  1090. Purpose:   Turbo Pascal does not have a arc cosine
  1091.  
  1092. Interface: function arccos( x : single) : single;
  1093.  
  1094. Errors:    Input value x must be within plus and minus 1.0
  1095.  
  1096. Remarks:   Uses arctan identity to solve for arccos
  1097.  
  1098. Example:   begin
  1099.                 x = rad_deg(arccos(cos(deg_rad(45.0))));  ( x = 45.0)
  1100.            end;
  1101.  
  1102.  
  1103. 4.3.3 Arc Tangent : four quadrants
  1104.  
  1105. Purpose:   Turbo Pascal does not have a arc tangent for all four 
  1106.            quadrants
  1107.  
  1108. Interface: function arctan2( x,y : single) : single;
  1109.  
  1110. Errors:    x and y values must not be zero
  1111.  
  1112. Remarks:   Uses Turbo Pascal atan to initially solve the arctan2 and 
  1113.            then the signs of the x and y terms dictate the quadrant
  1114.  
  1115. Example:   begin
  1116.             x := rad_deg(arctan2(tan(deg_rad(45.0),1.0))); ( x = 45.0)
  1117.            end;
  1118. 4.4 Hyperbolic Functions
  1119.  
  1120.  
  1121. 4.4.1 Hyperbolic sine
  1122.  
  1123. Purpose:   Turbo Pascal does not have a sinh function
  1124.  
  1125. Interface: function sinh( x : single) : single;
  1126.  
  1127. Errors:    None
  1128.  
  1129. Remarks:   None
  1130.  
  1131. Example:   begin
  1132.                 x = sinh(5.0)   (x = 74.210)
  1133.            end;
  1134.  
  1135.  
  1136. 4.4.2 Hyperbolic Cosine
  1137.  
  1138. Purpose:   Turbo Pascal does not have a cosh function
  1139.  
  1140. Interface: function cosh( x : single) : single;
  1141.  
  1142. Errors:    None
  1143.  
  1144. Remarks:   None
  1145.  
  1146. Example:   begin
  1147.                 x = cosh(5.0)   (x = 74.203)
  1148.            end;
  1149.  
  1150.  
  1151. 4.4.3 Hyperbolic Tangent
  1152.  
  1153. Purpose:   Turbo Pascal does not have a tanh function
  1154.  
  1155. Interface: function tanh( x : single) : single;
  1156.  
  1157. Errors:    None
  1158.  
  1159. Remarks:   None
  1160.  
  1161. Example:   begin
  1162.                 x = tanh(5.0)   (x = 1.0)
  1163.            end;
  1164. 4.5 Logarithmic Functions
  1165.  
  1166.  
  1167. 4.5.1 Power
  1168.  
  1169. Purpose:   To raise a real number to a real number (e.g., x**y)
  1170.  
  1171. Interface: function power(x,y: extended) : extended;
  1172.            Where x = number to be raised
  1173.                  y = raised value
  1174.  
  1175. Errors:    None
  1176.  
  1177. Remarks:   Uses standard logarithmic functions to perform function. 
  1178.            Extended precision is used to increase precision of 
  1179.            resulting number and increase the range of answer.
  1180.  
  1181. Example:   begin
  1182.                 z := power(2.0,3.0)  ( z = 8.0)
  1183.            end;
  1184.  
  1185.  
  1186. 4.5.2 Logarithms to base 10
  1187.  
  1188. Purpose:   Turbo Pascal does not provide a logarithms to base 10
  1189.  
  1190. Interface: function log10(x : single) : single;
  1191.  
  1192. Errors:    None
  1193.  
  1194. Remarks:   Uses standard logarithmic functions to perform function
  1195.  
  1196. Example:   begin
  1197.                 z := log10(1000)  ( z = 3.0)
  1198.            end;
  1199.  
  1200.  
  1201. 4.5.3 Logarithms to any base
  1202.  
  1203. Purpose:   Turbo Pascal does not provide a logarithms to any base
  1204.  
  1205. Interface: function logxy(x,y : single) : single;
  1206.            Where y = base
  1207.                  x = base 10 number to be converted to base y
  1208.            
  1209. Errors:    ln(y) must not be zero
  1210.  
  1211. Remarks:   Uses standard logarithmic functions to perform function
  1212.  
  1213. Example:   begin
  1214.                 z := logxy(1000.0,2.0)  ( z = 9.966)
  1215.            end;
  1216. 4.6 Other Useful Mathematical Functions
  1217.  
  1218.  
  1219. 4.6.1  Double Precision Product
  1220.  
  1221. Purpose:   It is very difficult in Turbo Pascal to guarantee that two 
  1222.            single precision number can be correctly be multiplied to 
  1223.            yield a double precision answer 
  1224.            
  1225. Interface: function dprod( x,y: extended) : extended;
  1226.  
  1227. Errors:    None
  1228.  
  1229. Remarks:   The two single precision numbers may be actual parameters 
  1230.            of the function call.  The formal arguments are converted 
  1231.            to extended precision, multiplied in extended precision and 
  1232.            then returned in extended precision.
  1233.  
  1234. Example:   var x,y : single;
  1235.                z   : extended;
  1236.            begin
  1237.                 z := dprod(1/3,1.0);  (z = 0.333333333333333333)
  1238.                                       (not z =0.333333343 etc)
  1239.            end; 
  1240.  
  1241.  
  1242. 4.6.2  Convert to Double Precision
  1243.  
  1244. Purpose:   It is very difficult in Turbo Pascal to convert a single  
  1245.            precision number to a double precision number. 
  1246.            
  1247. Interface: function dble( x extended) : extended;
  1248.  
  1249. Errors:    None
  1250.  
  1251. Remarks:   Not that there is no gain in precision of the original 
  1252.            number.
  1253.  
  1254. Example:   var x   : single;
  1255.                z   : extended;
  1256.            begin
  1257.                 z := dble(1/3);  
  1258.            end; 
  1259.  
  1260. 4.6.3  Root of an Equation : secant method
  1261.  
  1262. Purpose: Determine the root of an equation using secant method [3]
  1263.  
  1264. Interface: procedure secantmethod( xn, xn_1, fn, fn_1 : extended) : 
  1265.                                     extended;
  1266.            Where xn   = x(n)
  1267.                  nn_1 = x(n-1)
  1268.                  fn   = f(n) 
  1269.                  fn_1 = f(n_1)
  1270.  
  1271. Errors:    None
  1272.  
  1273. Remarks:   None 
  1274.  
  1275. Example:var 
  1276.            xn,xn_1,fxn,fxn_1 : extended; 
  1277.            guess : integer; 
  1278.  
  1279.            function formula (xn: extended) : extended; 
  1280.            const 
  1281.                 c2 : extended = 2.0; 
  1282.                 c3 : extended = 3.0; 
  1283.            begin 
  1284.                 formula := power(xn,c3) - c2*power(xn,c2) + xn - c3; 
  1285.            end; 
  1286.            begin 
  1287.                 xn   := 3.0; 
  1288.                 xn_1 := 4.0; 
  1289.                 repeat 
  1290.                      fxn   := formula(xn); 
  1291.                      fxn_1 := formula(xn_1); 
  1292.                      secantmethod(xn,xn_1,fxn,fxn_1); 
  1293.                 until abs(xn - xn_1) <= 1.0E-17; 
  1294.            converges to 2.17455941029297 etc.
  1295.            end.
  1296.  
  1297.  
  1298. 4.6.4  Factorial
  1299.  
  1300. Purpose:   Provide a n! function
  1301.  
  1302. Interface: function factorial( n: word) : single;
  1303.  
  1304. Errors:    None
  1305.  
  1306. Remarks:   None
  1307.  
  1308. Example:   begin
  1309.                 z := factorial(5);  (z = 120)
  1310.            end;
  1311. 5.0 References
  1312.  
  1313. [1] Mastering Turbo Pascal Version 5.5, Third Edition, Tom Swan, 
  1314. Hayden Books, 1990
  1315.  
  1316. [2] Pascal Programs : For Scientists and Engineers, Alan Miller, 
  1317. Sybex, 1981
  1318.  
  1319. [3] Numerical Mathematics and Computing, Ward Cheney and David 
  1320. Kincaid, Brooks/Cole Publishing, 1980
  1321.  
  1322. [4] SPSS/PC+ User's Manual, SPSS Inc.,  1986
  1323.  
  1324. [5] Digital Computer Methods in Engineering, S. Hovanessian and 
  1325. L.Pipes, McGraw-Hill, 1969
  1326.  
  1327. [6] Programming in Standard FORTRAN 77, A. Balfour and D. Marwick, 
  1328. North-Holland, 1979
  1329.  
  1330. [7] Introduction to Statistical Analysis, 
  1331.  
  1332. [8] Statistical Methods, H. Arkin and R. Colton, Barnes & Noble, 1970
  1333.  
  1334. [9] Data Reduction and Error Analysis for the Physical Sciences, P. 
  1335. Bevington, McGraw-Hill, 1969
  1336.  
  1337. [10] Calculus with Analytical Geometry, E. Purcell, 
  1338. Appleton-Century-Crofts, 1965
  1339.  
  1340. [11] Knuth, Volume 2 
  1341.  
  1342. [12] BYTE Magazine, March 1987
  1343.  
  1344. [13] Handbook of Statistical Methods in Meteorology, 
  1345.  
  1346. [14] Handbook of Mathematical Functions, Ambrowitz, 1982
  1347.