DECLARE SUB NormalArray (a!(), n%) DECLARE SUB NormalStats (a!(), n%, flag AS STRING, avg!, stddev!) DECLARE SUB LinearReg (x!(), y!(), n%, flag AS STRING, a!, b!, syx!, r!) DECLARE FUNCTION DotProduct! (n%, x!(), y!()) DECLARE FUNCTION Sum! (n%, x!()) '******************************************************************** '* This Program tests Module STATS for basic statistical Functions * '* ---------------------------------------------------------------- * '* SAMPLE RUN: * '* * '* Population mean and standard deviation: .0834715 1.000265 * '* * '* Sample mean and standard deviation: .0834715 1.005268 * '* * '* For full regression: * '* regression coefficients: -1.215311 2.048838 * '* standard error of estimate of y on x: 10.5399 * '* correlation coefficient: .9841678 * '* * '* For regression forced through (0,0): * '* regression coefficients: 0 2.024772 * '* standard error of estimate of y on x: 10.56329 * '* correlation coefficient: .9840971 * '* * '* ---------------------------------------------------------------- * '* Ref.: "Problem Solving with Fortran 90 By David R.Brooks, * '* Springer-Verlag New York, 1997". * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '******************************************************************** 'Program Tstats DEFINT I-N n = 100 DIM x(n), y(n) ' Main variables ' avg, stddev average, standard deviation ' a, b, for linear regression y=ax+b ' syx standard error of estimate of y on x ' corr correlation coefficient ' n # of points in arrays x, y ' Generate random array x of size n NormalArray x(), n ' Test basic statistics on array x NormalStats x(), n, "p", avg, stddev CLS PRINT PRINT " Population mean and standard deviation: "; avg; " "; stddev NormalStats x(), n, "s", avg, stddev PRINT PRINT " Sample mean and standard deviation: "; avg; " "; stddev PRINT ' Test linear regression ' Generate random array y of size n NormalArray y(), n ' Create a linear relationship with "noise" FOR i = 1 TO n x(i) = 1! * i y(i) = 2! * i + 10! * y(i) NEXT i ' Set a <> zero for full regression analysis a = 1! LinearReg x(), y(), n, "s", a, b, syx, corr PRINT " For full regression:" PRINT " regression coefficients: "; a; " "; b PRINT " standard error of estimate of y on x: "; syx PRINT " correlation coefficient: "; corr PRINT ' Set a = zero for full regression forced through (0,0) a = 0! LinearReg x(), y(), n, "s", a, b, syx, corr PRINT " For regression forced through (0,0):" PRINT " regression coefficients: "; a; " "; b PRINT " standard error of estimate of y on x: "; syx PRINT " correlation coefficient: "; corr PRINT END 'of main program FUNCTION DotProduct (n, x(), y()) temp = 0! FOR i = 1 TO n temp = temp + x(i) * y(i) NEXT i DotProduct = temp END FUNCTION 'DotProduct SUB LinearReg (x(), y(), n, flag AS STRING, a, b, syx, r) '---------------------------------------------------------- ' For data to be represented by y=ax+b, calculates linear ' regression coefficients, sample standard error of y on x, ' and sample correlation coefficients. Sets r=0 if an error ' exists. If the intercept coefficient a is set to 0 on ' input, the regression is forced through (0,0). '---------------------------------------------------------- ' avg, stddev: Single ' sumx,sumy,sumxy,sumxx,sumyy,temp: Single sumx = Sum(n, x()) sumy = Sum(n, y()) sumxy = DotProduct(n, x(), y()) sumxx = DotProduct(n, x(), x()) sumyy = DotProduct(n, y(), y()) IF a <> 0! THEN 'calculate full expression temp = n * sumxx - sumx ^ 2 a = (sumy * sumxx - sumx * sumxy) / temp b = (n * sumxy - sumx * sumy) / temp syx = SQR((sumyy - a * sumy - b * sumxy) / n) ELSE 'just calculate slope b = sumy / sumx syx = SQR((sumyy - 2! * b * sumxy + b * b * sumxx) / n) END IF IF flag = "s" OR flag = "S" THEN syx = syx * SQR((1! * n) / (1! * (n - 2))) END IF ' Use NormalStats to get standard deviation of y NormalStats y(), n, flag, avg, stddev IF stddev > 0! THEN temp = 1! - (syx / stddev) ^ 2 IF temp > 0! THEN r = SQR(temp) ELSE 'an error exists r = 0! PRINT " From LinearReg: error in temp "; temp END IF ELSE 'an error exists r = 0! END IF END SUB 'LinearReg SUB NormalArray (a(), n) '----------------------------------------------------------- ' Generates an array of normal random numbers from pairs of ' uniform random numbers in range [0,1]. '----------------------------------------------------------- pi = 3.141593 'fills array with uniform random FOR i = 1 TO n a(i) = RND NEXT i FOR i = 1 TO n STEP 2 u1 = a(i) u2 = a(i + 1) IF u1 = 0! THEN u1 = 1E-12 'u must not be zero IF u2 = 0! THEN u2 = 1E-12 a(i) = SQR(-2! * LOG(u1)) * COS(2! * pi * u2) a(i + 1) = SQR(-2! * LOG(u2)) * SIN(2! * pi * u2) NEXT i IF (n MOD 2) <> 0! THEN 'there is one extra element IF a(n) = 0! THEN a(n) = 1E-12 a(n) = SQR(-2! * LOG(a(n))) * SIN(2! * pi * a(n)) END IF END SUB 'NormalArray SUB NormalStats (a(), n, flag AS STRING, avg, stddev) '--------------------------------------------------------- ' Basic description statistics for normally distributed ' data. Calculates either sample or population statistics. ' Sets stddev to -1 if error condition is detected. '--------------------------------------------------------- ' sum0, sumsq, variance: Single sum0 = Sum(n, a()) sumsq = DotProduct(n, a(), a()) IF flag = "p" OR flag = "P" THEN variance = (sumsq - sum0 ^ 2 / n) / n ELSEIF flag = "s" OR flag = "S" THEN variance = (sumsq - sum0 ^ 2 / (n - 1)) / (n - 1) ELSE PRINT " From NormalStats: Flag Error,

assumed." variance = (sumsq - sum0 ^ 2 / n) / n END IF IF variance < 0! THEN 'an error exists PRINT " From NormalStats: negative variance "; variance stddev = -1! ELSE stddev = SQR(variance) END IF avg = sum0 / n END SUB 'NormalStats FUNCTION Sum (n, x()) ' Returns sum of all elements of array x temp = 0! FOR i = 1 TO n temp = temp + x(i) NEXT i Sum = temp END FUNCTION 'Sum