21 #include "./base/base_uses.f90"
26 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'statistical_methods'
27 LOGICAL,
PARAMETER :: debug_this_module = .false.
47 REAL(kind=
dp),
DIMENSION(:),
POINTER :: ix
48 INTEGER,
INTENT(IN) :: n
49 REAL(kind=
dp),
INTENT(OUT) :: w, pw
51 REAL(kind=
dp),
PARAMETER :: c1(6) = (/0.000000_dp, 0.221157_dp, -0.147981_dp, -2.071190_dp, &
52 4.434685_dp, -2.706056_dp/), c2(6) = (/0.000000_dp, 0.042981_dp, -0.293762_dp, &
53 -1.752461_dp, 5.682633_dp, -3.582633_dp/), &
54 c3(4) = (/0.544000_dp, -0.399780_dp, 0.025054_dp, -0.6714e-3_dp/), &
55 c4(4) = (/1.3822_dp, -0.77857_dp, 0.062767_dp, -0.0020322_dp/), &
56 c5(4) = (/-1.5861_dp, -0.31082_dp, -0.083751_dp, 0.0038915_dp/), &
57 c6(3) = (/-0.4803_dp, -0.082676_dp, 0.0030302_dp/), g(2) = (/-2.273_dp, 0.459_dp/), &
58 one = 1.0_dp, pi6 = 1.909859_dp, qtr = 0.25_dp, small = epsilon(0.0_dp), &
60 REAL(kind=
dp),
PARAMETER :: stqr = 1.047198_dp, th = 0.375_dp, two = 2.0_dp, zero = 0.0_dp
62 INTEGER :: i, i1, j, n2, output_unit
63 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: itmp
65 REAL(kind=
dp) :: a1, a2, an, an25, asa,
fac,
gamma, m, &
66 range, rsn, s, sa, sax, ssa, ssassx, &
67 ssumm2, ssx, summ2, sx, w1, xi, xsx, &
69 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: a, x
74 IF (n < 3 .OR. n > 5000)
THEN
75 IF (output_unit > 0)
WRITE (output_unit,
'(A)') &
76 "Shapiro Wilk test: number of points less than 3 or greated than 5000."
77 IF (output_unit > 0)
WRITE (output_unit,
'(A)') &
78 "Not able to perform the test!"
81 IF (mod(n, 2) == 0)
THEN
93 IF (range < small) failure = .true.
94 IF (failure .AND. (output_unit > 0))
THEN
95 WRITE (output_unit,
'(A)')
"Shapiro Wilk test: two data points are numerically identical."
96 WRITE (output_unit,
'(A)')
"Not able to perform the test!"
108 CALL ppnd7((i - th)/an25, a(i))
109 summ2 = summ2 + a(i)**2
114 a1 = poly(c1, 6, rsn) - a(1)/ssumm2
118 a2 = -a(2)/ssumm2 + poly(c2, 6, rsn)
119 fac = sqrt((summ2 - two*a(1)**2 - two*a(2)**2)/(one - two*a1**2 - two*a2**2))
124 fac = sqrt((summ2 - two*a(1)**2)/(one - two*a1**2))
139 IF (i /= j) sa = sa + sign(1, i - j)*a(min(i, j))
153 asa = sign(1, i - j)*a(min(i, j)) - sa
157 xsx = x(i)/range - sx
165 ssassx = sqrt(ssa*ssx)
166 w1 = (ssassx - sax)*(ssassx + sax)/(ssa*ssx)
170 pw = pi6*(asin(sqrt(w)) - stqr)
177 gamma = poly(g, 2, an)
183 s = exp(poly(c4, 4, an))
184 pw = alnorm((y - m)/s, .true.)
188 s = exp(poly(c6, 3, xx))
189 pw = alnorm((y - m)/s, .true.)
208 SUBROUTINE ppnd7(p, normal_dev)
209 REAL(kind=
dp),
INTENT(IN) :: p
210 REAL(kind=
dp),
INTENT(OUT) :: normal_dev
212 REAL(kind=
dp),
PARAMETER :: a0 = 3.3871327179e+00_dp, a1 = 5.0434271938e+01_dp, &
213 a2 = 1.5929113202e+02_dp, a3 = 5.9109374720e+01_dp, b1 = 1.7895169469e+01_dp, &
214 b2 = 7.8757757664e+01_dp, b3 = 6.7187563600e+01_dp, c0 = 1.4234372777e+00_dp, &
215 c1 = 2.7568153900e+00_dp, c2 = 1.3067284816e+00_dp, c3 = 1.7023821103e-01_dp, &
216 const1 = 0.180625_dp, const2 = 1.6_dp, d1 = 7.3700164250e-01_dp, &
217 d2 = 1.2021132975e-01_dp, e0 = 6.6579051150e+00_dp, e1 = 3.0812263860e+00_dp, &
218 e2 = 4.2868294337e-01_dp, e3 = 1.7337203997e-02_dp, f1 = 2.4197894225e-01_dp, &
219 f2 = 1.2258202635e-02_dp, half = 0.5_dp, one = 1.0_dp
220 REAL(kind=
dp),
PARAMETER :: split1 = 0.425_dp, split2 = 5.0_dp, zero = 0.0_dp
222 REAL(kind=
dp) :: q, r
225 IF (abs(q) <= split1)
THEN
227 normal_dev = q*(((a3*r + a2)*r + a1)*r + a0)/ &
228 (((b3*r + b2)*r + b1)*r + one)
241 IF (r <= split2)
THEN
243 normal_dev = (((c3*r + c2)*r + c1)*r + c0)/((d2*r + d1)*r + one)
246 normal_dev = (((e3*r + e2)*r + e1)*r + e0)/((f2*r + f1)*r + one)
248 IF (q < zero) normal_dev = -normal_dev
265 FUNCTION alnorm(x, upper)
RESULT(fn_val)
266 REAL(kind=
dp),
INTENT(IN) :: x
267 LOGICAL,
INTENT(IN) :: upper
268 REAL(kind=
dp) :: fn_val
270 REAL(kind=
dp),
PARAMETER :: a1 = 5.75885480458_dp, a2 = 2.62433121679_dp, &
271 a3 = 5.92885724438_dp, b1 = -29.8213557807_dp, b2 = 48.6959930692_dp, c1 = -3.8052e-8_dp, &
272 c2 = 3.98064794e-4_dp, c3 = -0.151679116635_dp, c4 = 4.8385912808_dp, &
273 c5 = 0.742380924027_dp, c6 = 3.99019417011_dp, con = 1.28_dp, d1 = 1.00000615302_dp, &
274 d2 = 1.98615381364_dp, d3 = 5.29330324926_dp, d4 = -15.1508972451_dp, &
275 d5 = 30.789933034_dp, half = 0.5_dp, ltone = 7.0_dp, one = 1.0_dp, p = 0.398942280444_dp, &
276 q = 0.39990348504_dp, r = 0.398942280385_dp, utzero = 18.66_dp, zero = 0.0_dp
279 REAL(kind=
dp) :: y, z
287 IF (.NOT. (z <= ltone .OR.
up .AND. z <= utzero))
THEN
289 IF (.NOT.
up) fn_val = one - fn_val
294 fn_val = r*exp(-y)/(z + c1 + d1/(z + c2 + d2/(z + c3 + d3/(z + c4 + d4/(z + c5 + d5/(z + c6))))))
296 fn_val = half - z*(p - q*y/(y + a1 + b1/(y + a2 + b2/(y + a3))))
298 IF (.NOT.
up) fn_val = one - fn_val
314 FUNCTION poly(c, nord, x)
RESULT(fn_val)
316 REAL(kind=
dp),
INTENT(IN) :: c(:)
317 INTEGER,
INTENT(IN) :: nord
318 REAL(kind=
dp),
INTENT(IN) :: x
319 REAL(kind=
dp) :: fn_val
325 IF (nord == 1)
RETURN
355 SUBROUTINE k_test(xdata, istart, n, tau, z, prob)
356 REAL(kind=
dp),
DIMENSION(:),
POINTER :: xdata
357 INTEGER,
INTENT(IN) :: istart, n
358 REAL(kind=
dp) :: tau, z, prob
360 INTEGER :: is, j, k, nt
361 REAL(kind=
dp) :: a1, var
368 a1 = xdata(j) - xdata(k)
369 IF (a1 .GT. 0.0_dp)
THEN
371 ELSE IF (a1 .LT. 0.0_dp)
THEN
376 tau = real(is, kind=
dp)
377 var = real(nt, kind=
dp)*real(nt - 1, kind=
dp)*real(2*nt + 5, kind=
dp)/18.0_dp
379 prob = erf(abs(z)/sqrt(2.0_dp))
398 REAL(kind=
dp),
DIMENSION(:),
POINTER :: xdata
399 INTEGER,
INTENT(IN) :: n
400 REAL(kind=
dp) :: r, u, prob
403 REAL(kind=
dp) :: q, s, var, x
411 q = q + (xdata(i + 1) - xdata(i))**2
414 x = x/real(n, kind=
dp)
416 s = s + (xdata(i) - x)**2
418 s = s/real(n - 1, kind=
dp)
419 q = q/real(2*(n - 1), kind=
dp)
421 var = sqrt(1.0_dp/real(n + 1, kind=
dp)*(1.0_dp + 1.0_dp/real(n - 1, kind=
dp)))
423 prob = erf(abs(u)/sqrt(2.0_dp))
440 SUBROUTINE tests(xdata, globenv)
441 REAL(kind=
dp),
DIMENSION(:),
POINTER :: xdata
442 TYPE(global_environment_type),
POINTER :: globenv
445 REAL(kind=
dp) :: prob, pw, r, tau, u, w, z
446 REAL(kind=
dp),
DIMENSION(:),
POINTER :: ydata
448 IF (debug_this_module)
THEN
453 xdata(i) = 5.0_dp - real(i, kind=
dp)/2.0_dp + 0.1*globenv%gaussian_rng_stream%next()
454 WRITE (3, *) xdata(i)
457 xdata(i) = 0.1*globenv%gaussian_rng_stream%next()
462 CALL k_test(xdata, i, n, tau, z, prob)
463 IF (prob <= 0.2_dp)
EXIT
465 WRITE (*, *)
"Mann-Kendall test", i
469 ALLOCATE (ydata(n - i + 1))
471 CALL sw_test(ydata, n - i + 1, w, pw)
472 CALL vn_test(ydata, n - i + 1, r, u, prob)
473 WRITE (*, *)
"Shapiro Wilks test", i, w, pw
474 WRITE (*, *)
"Von Neu", i, r, u, prob
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
static GRID_HOST_DEVICE orbital up(const int i, const orbital a)
Increase i'th component of given orbital angular momentum.
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Define type storing the global information of a run. Keep the amount of stored data small....
Defines the basic variable types.
integer, parameter, public dp
Methods to perform on the fly statistical analysis of data -) Schiferl and Wallace,...
subroutine, public vn_test(xdata, n, r, u, prob)
Von Neumann test for serial correlation.
integer, parameter, public min_sample_size
subroutine, public k_test(xdata, istart, n, tau, z, prob)
Kandall's test for correlation.
subroutine, public sw_test(ix, n, w, pw)
Shapiro - Wilk's test or W-statistic to test normality of a distribution R94 APPL....
All kind of helpful little routines.