!****************************************** !* Module for basic statistical Functions * !****************************************** Module Stats Contains !--------------------------------------------------------- ! Basic description statistics for normally distributed ! data. Calculates either sample or population statistics. ! Sets std_dev to -1 if error condition is detected. !--------------------------------------------------------- Subroutine NormalStats(a,n,flag,avg,std_dev) Implicit None Integer, Intent(IN) :: n Real, Intent(IN) :: a(n) Character, Intent(IN) :: flag ! 'P' or 'S' statistics Real, Intent(OUT) :: avg, std_dev Real sum_, sum_sq, variance sum_ = SUM(a) sum_sq = DOT_PRODUCT(a,a) Select Case(flag) Case ('p','P') variance = (sum_sq-sum_**2/n)/n Case ('s','S') variance = (sum_sq-sum_**2/n)/(n-1) Case DEFAULT print *,' From NormalStats: Flag Error,

assumed.' variance = (sum_sq-sum_**2/n)/n End Select If (variance < 0.) Then !an error exists print *,' From NormalStats: negative variance ',variance std_dev = -1 Else std_dev = SQRT(variance) End If avg = sum_/n End Subroutine NormalStats !---------------------------------------------------------- ! 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). !---------------------------------------------------------- Subroutine LinearReg(x,y,n,flag,a,b,s_yx,r) Implicit None Integer, Intent(IN) :: n Real, Intent(IN) :: x(n), y(n) Character, Intent(IN) :: flag ! 'P' or 'S' statistics Real, Intent(OUT) :: b, s_yx, r Real, Intent(INOUT) :: a Real avg, std_dev Real sum_x,sum_y,sum_xy,sum_xx,sum_yy,temp sum_x = SUM(x) sum_y = SUM(y) sum_xy = DOT_PRODUCT(x,y) sum_xx = DOT_PRODUCT(x,x) sum_yy = DOT_PRODUCT(y,y) If (a /= 0.) Then !calculate full expression temp = n*sum_xx - sum_x**2 a = (sum_y*sum_xx - sum_x*sum_xy)/temp b = (n*sum_xy - sum_x*sum_y)/temp s_yx = SQRT((sum_yy - a*sum_y - b*sum_xy)/n) Else !just calculate slope b = sum_y/sum_x s_yx = SQRT((sum_yy - 2.*b*sum_xy + b*b*sum_xx)/n) End If Select Case(flag) Case ('p','P') Case ('s','S') s_yx = s_yx * SQRT(Real(n)/Real(n-2)) Case DEFAULT print *,' From LinearReg: Flag error,

assumed.' End Select ! Use NormalStats to get standard deviation of y Call NormalStats(y,n,flag,avg,std_dev) If (std_dev > 0.) Then temp = 1. - (s_yx/std_dev)**2 If (temp > 0.) Then r = SQRT(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 Subroutine LinearReg !----------------------------------------------------------- ! Generates an array of normal random numbers from pairs of ! of uniform random numbers in range [0,1]. !----------------------------------------------------------- Subroutine NormalArray(a,n) Implicit None Integer, Intent(IN) :: n Real, Intent(OUT) :: a(n) Integer i Integer Count(1) !for random number generator Real pi, u1,u2 Parameter(pi=3.1415926535) Call System_Clock(Count(1)) Call Random_Seed(Put=Count) Call Random_Number(a) !fills array with uniform random Do i=1,n,2 u1 = a(i) u2 = a(i+1) If (u1==0.) u1 = 1e-15 !u must not be zero If (u2==0.) u2 = 1e-15 a(i) = SQRT(-2.0*LOG(u1))*COS(2.0*pi*u2) a(i+1) = SQRT(-2.0*LOG(u2))*SIN(2.0*pi*u2) End Do If (MOD(n,2) /= 0) Then !there is one extra element If (a(n) == 0.) a(n) = 1e-15 a(n) = SQRT(-2.0*LOG(a(n)))*SIN(2.0*pi*a(n)) End If End Subroutine NormalArray End Module Stats !end of file Stats.f90