chemshell_qchem / blob
1 PROGRAM protTool
2 !==============================================================================
3 ! Reads a data set, groups the data into bins, calculates bin means and standard
4 ! deviations, and performs a series of statistical tests on the coarse-grained
5 ! data. Tests are for trend, normality, and serial correlation.
6 ! Procedure follows
7 ! Schiferl, Wallace, J. Chem. Phys. 1985, 83, 5203.
8 !
9 ! Hans Martin Senn, MPI for Coal Research, Muelheim an der Ruhr, 2003.
10 !------------------------------------------------------------------------------
11 ! USAGE: protTool <file> <column #> <skip> <bin width> [<result file>]
12 ! Input data is assumed to be arranged in columns; lines starting with # are skipped.
13 !------------------------------------------------------------------------------
14 ! DEPENDENCIES
15 ! We use the following routines from the Compaq cxml library:
16 ! - dnrsq: Sum of squares of the components of a vector.
17 ! - dsortq: Sorts a vector.
18 ! We use the following routines from the fortran(3f) library:
19 ! - derfc: Error function
20 !==============================================================================
21 IMPLICIT NONE
23 ! File units
24 INTEGER(4), PARAMETER :: protFile=10, resFile=20
26 INTEGER(4), PARAMETER :: Ntests = 4
28 INTEGER(4) :: col, skip, width, Nlines, Nval, Nbins, excess
29 INTEGER(4) :: i
30 REAL(8) :: meanSD(2), CI
31 REAL(8), ALLOCATABLE :: val(:), stat(:,:)
32 CHARACTER(12) :: res(Ntests) = 'NOT PASSED'
33 CHARACTER(40) :: test(Ntests)
34 LOGICAL(4) :: lres, Tpass(Ntests) = .FALSE.
35 ! ------------------------------------------------------------------------------
37 ! Read command line input, open files
38 CALL input(protFile, resFile, lres, col, skip, width, Nlines)
40 ! Read data
41 ALLOCATE(val(Nlines))
42 CALL read_data(protFile, col, val, Nval, Nlines)
43 CLOSE(protFile)
45 ! Arrange data into bins
46 Nbins = (Nval - skip) / width
47 excess = (Nval - skip) - Nbins*width
48 WRITE(*,*)
49 WRITE(*,*) 'No. of bins: ', Nbins
50 WRITE(*,*) 'No. of data points discarded to allow regular binning: ', excess
51 skip = skip + excess
52 WRITE(*,*) 'Adjusted no. of skipped values: ', skip
53 WRITE(*,*) '-----------------------------------------------------------------'
55 ! Mean and SD in each bin
56 ALLOCATE(stat(Nbins,2))
57 WRITE(*,*)
58 WRITE(*,*) 'Calculating mean and standard deviation of each bin...'
59 CALL mean_sd(val, stat, skip, width, Nval, Nbins) ! Means and SD's are in stat
61 IF (lres) THEN
62 WRITE(resFile,'(I10, 2G20.8)') &
63 & (skip + i*width - width/2, stat(i, 1), stat(i, 2), i = 1,Nbins)
64 CLOSE(resFile)
65 END IF
67 ! Coarse-grained mean and SD (i.e., mean and SD of bin means)
68 CALL mean_sd(stat(:,1), meanSD, 0, Nbins, Nbins, 1)
70 ! Confidence interval at 95%, for N > 24 (Student t fractile approx. as 2)
71 CI = 2.d0*meanSD(2)/SQRT(DBLE(Nbins))
73 ! Test for trend in means and variances
74 WRITE(*,*)
75 WRITE(*,*) 'Testing for TREND IN BIN MEANS...'
76 CALL Mann_Kendall(stat(:,1), Nbins, Tpass(1))
78 WRITE(*,*)
79 WRITE(*,*) 'Testing for TREND IN BIN STANDARD DEVIATIONS...'
80 CALL Mann_Kendall(stat(:,2), Nbins, Tpass(2))
82 ! Test for normality
83 WRITE(*,*)
84 WRITE(*,*) 'Testing for NORMALITY OF BIN MEANS...'
85 CALL Shapiro_Wilk(stat(:,1), Nbins, Tpass(3))
87 ! Test for correlation
88 WRITE(*,*)
89 WRITE(*,*) 'Testing for POSITIVE CORRELATION OF BIN MEANS...'
90 CALL Neumann(stat(:,1), meanSD(2)**2, Nbins, Tpass(4))
92 ! The final verdict
93 test(1) = ' Lack of trend in bin means'
94 test(2) = ' Lack of trend in bin SD''s'
95 test(3) = ' Normality of bin means'
96 test(4) = ' Lack of correlation in bin means'
98 WHERE (Tpass) res = 'passed'
100 WRITE(*,*)
101 WRITE(*,*) '-----------------------------------------------------------------'
102 WRITE(*,*) '---------------------- SUMMARY ---------------------------'
103 WRITE(*,*) '-----------------------------------------------------------------'
104 WRITE(*,*)
105 WRITE(*,*) 'Number of observations: ', Nbins
106 WRITE(*,'((A),ES12.4)') ' Sample mean: ', meanSD(1)
107 WRITE(*,'((A),ES12.4)') ' Confidence limits at 95%: +/-', CI
108 WRITE(*,*) '(NB: Valid for N > 24)'
109 WRITE(*,*)
110 WRITE(*,'(((A), (A)))') (test(i), TRIM(res(i)), i=1,Ntests)
112 WRITE(*,*) '-----------------------------------------------------------------'
115 DEALLOCATE(val)
116 DEALLOCATE(stat)
117 END PROGRAM protTool
120 !==============================================================================
121 SUBROUTINE input(protFile, resFile, lres, col, skip, width, Nlines)
122 !==============================================================================
123 IMPLICIT NONE
125 INTEGER(4),INTENT(in) :: protFile, resFile
126 INTEGER(4),INTENT(out) :: skip, width, col, Nlines
127 LOGICAL(4),INTENT(out) :: lres
129 INTEGER(4) :: iargc, Nargs
130 INTEGER(4) :: ios = 0
131 CHARACTER(64) :: arg, protName, resName
132 CHARACTER(128) :: line
134 ! Compaq fortran(3f) routines
135 EXTERNAL getarg, iargc
136 ! ------------------------------------------------------------------------------
137 ! Parse command line
138 ! ------------------------------------------------------------------------------
139 Nargs = iargc()
140 IF (Nargs < 4) THEN
141 WRITE(*,*)
142 WRITE(*,*) 'USAGE:'
143 WRITE(*,*) 'protTool <prot file> <column #> <skip> <bin width> [<result file>]'
144 WRITE(*,*)
145 STOP
146 END IF
148 CALL getarg(1, protName)
150 CALL getarg(2, arg)
151 arg = TRIM(arg)
152 READ(arg,*) col
154 CALL getarg(3, arg)
155 arg = TRIM(arg)
156 READ(arg,*) skip
158 CALL getarg(4, arg)
159 arg = TRIM(arg)
160 READ(arg,*) width
162 IF (Nargs > 4) THEN
163 CALL getarg(5, arg)
164 arg = TRIM(arg)
165 READ(arg,*) resName
166 lres = .TRUE.
167 ELSE
168 resName = ''
169 lres = .FALSE.
170 END IF
172 ! Report input parameters
173 WRITE(*,*)
174 WRITE(*,*) '================================================================='
175 WRITE(*,*) '=================== protTool ==================='
176 WRITE(*,*) '================================================================='
177 WRITE(*,*)
178 WRITE(*,*) '-----------------------------------------------------------------'
179 WRITE(*,*) '------------------- Input Parameters -------------------'
180 WRITE(*,*) '-----------------------------------------------------------------'
181 WRITE(*,*) 'Data file: ', protName
182 WRITE(*,*) 'Results written to file: ', resName
183 WRITE(*,*) 'Data taken from column no.: ', col
184 WRITE(*,*) 'No. of steps skipped at the beginning: ', skip
185 WRITE(*,*) 'Bin width for coarse-graining: ', width
186 WRITE(*,*) '-----------------------------------------------------------------'
187 WRITE(*,*)
189 ! ------------------------------------------------------------------------------
190 ! Open result file if required
191 ! ------------------------------------------------------------------------------
192 IF (lres) THEN
193 OPEN(UNIT=resFile, FILE=TRIM(resName), STATUS='NEW', FORM='FORMATTED',&
194 & POSITION='REWIND', ACTION='WRITE')
195 END IF
197 ! ------------------------------------------------------------------------------
198 ! Open protocol for reading and estimate max. no. of values
199 ! ------------------------------------------------------------------------------
200 OPEN(UNIT=protFile, FILE=TRIM(protName), STATUS='OLD', FORM='FORMATTED',&
201 & POSITION='REWIND', ACTION='READ')
203 ! Count no. of lines in file
204 ! N.B.: Blank lines are not read (and not counted).
205 Nlines = 0
206 DO WHILE (ios == 0)
207 READ(protFile, *, IOSTAT=ios) line
208 SELECT CASE(ios)
209 CASE (:-1)
210 EXIT
211 CASE (1:)
212 STOP 'ERROR while reading prot file.'
213 END SELECT
214 Nlines = Nlines + 1
215 END DO
217 WRITE(*,*) '-----------------------------------------------------------------'
218 WRITE(*,*) '------------------- File & Binning Info -------------------'
219 WRITE(*,*) '-----------------------------------------------------------------'
220 WRITE(*,*) 'Non-empty lines in file: ', Nlines
222 RETURN
223 END SUBROUTINE input
226 !==============================================================================
227 SUBROUTINE read_data(protFile, col, val, Nval, Nlines)
228 !==============================================================================
229 IMPLICIT NONE
231 INTEGER(4), INTENT(in) :: protFile, col, Nlines
232 INTEGER(4), INTENT(out) :: Nval
233 REAL(8), INTENT(out) :: val(Nlines)
235 INTEGER(4) :: i
236 INTEGER(4) :: ios = 0
237 CHARACTER(64) :: record(col), ch64
238 ! ------------------------------------------------------------------------------
239 ! Read data from prot file
240 ! ------------------------------------------------------------------------------
241 REWIND(protFile)
242 Nval = 0
243 DO WHILE (ios == 0)
244 READ(protFile, *, IOSTAT=ios) (record(i), i=1,col)
245 SELECT CASE(ios)
246 CASE (:-1)
247 EXIT
248 CASE (1:)
249 STOP 'ERROR while reading prot file.'
250 END SELECT
251 ch64 = ADJUSTL(record(1))
252 IF (ch64(1:1).EQ.'#') THEN
253 CYCLE
254 ELSE
255 Nval = Nval + 1
256 READ(record(col), '(G64.12)') val(Nval)
257 END IF
258 END DO
260 WRITE(*,*) 'No. of values in data set: ', Nval
262 RETURN
263 END SUBROUTINE read_data
266 !==============================================================================
267 SUBROUTINE mean_sd(val, stat, skip, width, Nval, Nbins)
268 !==============================================================================
269 IMPLICIT NONE
271 INTEGER(4), INTENT(in) :: skip, width, Nval, Nbins
272 REAL(8), INTENT(in) :: val(Nval)
273 REAL(8), INTENT(out) :: stat(Nbins,2)
275 INTEGER(4) :: ibin
276 REAL(8) :: bin(width)
278 ! cxml function
279 REAL(8) :: dnrsq ! sum of squares
280 ! ------------------------------------------------------------------------------
281 ! Calculate mean and SD of each bin
282 ! ------------------------------------------------------------------------------
283 DO ibin=1,Nbins
284 bin = val(ibin*width + skip - width + 1:ibin*width + skip)
285 stat(ibin,1) = SUM(bin) / DBLE(width) ! mean
286 stat(ibin,2) = dnrsq(width, bin - stat(ibin,1), 1)
287 END DO
288 IF (width == 1) THEN
289 ! Catch pathological case of bin width = 1
290 stat(:,2) = 0.d0
291 ELSE
292 stat(:,2) = SQRT(stat(:,2) / DBLE(width - 1)) ! SD = sqrt(var)
293 END IF
295 RETURN
296 END SUBROUTINE mean_sd
299 !==============================================================================
300 SUBROUTINE Mann_Kendall(sample, Nobs, Tpass)
301 !==============================================================================
302 IMPLICIT NONE
304 INTEGER(4), INTENT(in) :: Nobs
305 REAL(8), INTENT(in) :: sample(Nobs)
306 LOGICAL(4), INTENT(out) :: Tpass
308 INTEGER(4) :: i, j, Inc, S
309 REAL(8) :: Ndble, SDS, u
311 ! Quantiles of standard normal distribution
312 REAL(8), PARAMETER :: u90 = 1.2816d0, &
313 & u95 = 1.6449d0, &
314 & u975 = 1.96d0, &
315 & u99 = 2.3263d0, &
316 & u995 = 2.5758d0
317 ! ------------------------------------------------------------------------------
318 ! Mann-Kendall test for trend
319 ! ------------------------------------------------------------------------------
320 Tpass = .FALSE.
321 WRITE(*,*) '-----------------------------------------------------------------'
322 WRITE(*,*) '-------------- Mann-Kendall test for trend --------------'
323 WRITE(*,*) '-----------------------------------------------------------------'
324 WRITE(*,*)
325 IF (Nobs < 20) THEN
326 WRITE(*,*) 'WARNING: The number of observations is <20.'
327 WRITE(*,*) ' The approximation of normal distribution for S is inaccurate!'
328 WRITE(*,*)
329 END IF
330 WRITE(*,*) 'Null hypothesis H0: There is no monotonic increase/decrease'
331 WRITE(*,*) ' with time in the sample (no trend).'
332 WRITE(*,*)
333 WRITE(*,*) 'Alternate hypothesis H1: There is a trend.'
334 WRITE(*,*)
335 WRITE(*,*) 'The test statistic S is approx. normally distributed with mean = 0'
336 WRITE(*,*) 'and standard deviation SD(S).'
337 WRITE(*,*) 'Assumption: No two observations are equal.'
338 WRITE(*,*)
340 ! Calculate no. of distinct increasing pairs of observations
341 ! (x_i, x_j) where x_j > x_i and j > i
342 Inc = 0
343 DO i = 1, Nobs - 1
344 DO j = i + 1, Nobs
345 IF (sample(j) > sample(i)) Inc = Inc + 1
346 END DO
347 END DO
349 ! Calculate the test statistic S
350 S = 2*Inc - Nobs*(Nobs - 1)/2 ! The second term is "n over 2"
352 ! SD of S
353 Ndble = DBLE(Nobs) ! N**3 needed below may well be too large for integer(4)
354 SDS = SQRT(Ndble*(Ndble - 1.d0)*(2.d0*Ndble + 5.d0) / 18.d0)
356 ! Standard normal deviate of S
357 u = DBLE(S)/SDS
359 WRITE(*,*) 'Test statistic S: ', S
360 WRITE(*,'((A),F10.3)') ' Standard deviation SD(S): ', SDS
361 WRITE(*,'((A),F10.3)') ' Standard normal deviate u = S/SD(S): ', u
362 WRITE(*,*) '(S lies u SD''s from the expected mean 0)'
363 WRITE(*,*)
365 WRITE(*,*) 'At significance level alpha = 0.05 (two-tailed test)'
366 IF (ABS(u) > u975) THEN
367 WRITE(*,*) '*** H0 IS REJECTED (there is a trend) ***'
368 WRITE(*,*) 'Recommendation: Increase start time (no. of skipped values).'
369 ELSE
370 WRITE(*,*) '*** H0 IS NOT REJECTED (no trend) ***'
371 Tpass = .TRUE.
372 END IF
373 WRITE(*,*)
375 WRITE(*,*) 'Quantiles of standard normal distribution:'
376 WRITE(*,*) 'alpha u(1-alpha) u(1-alpha/2)'
377 WRITE(*,'(F6.2,2F10.4)') 0.1, u90, u95
378 WRITE(*,'(F6.2,2F10.4)') 0.05, u95, u975
379 WRITE(*,'(F6.2,2F10.4)') 0.01, u99, u995
380 WRITE(*,*) '-----------------------------------------------------------------'
382 RETURN
383 END SUBROUTINE Mann_Kendall
386 !==============================================================================
387 SUBROUTINE Shapiro_Wilk(sample, Nobs, Tpass)
388 !==============================================================================
389 IMPLICIT NONE
391 INTEGER(4), INTENT(in) :: Nobs
392 REAL(8), INTENT(in) :: sample(Nobs)
394 LOGICAL(4), INTENT(out) :: Tpass
396 REAL(8), PARAMETER :: small = 1.d-24
398 REAL(8), DIMENSION(6), PARAMETER :: &
399 & C1 = (/0.0d0, 0.221157d0, -0.147981d0, -0.207119d1, 0.4434685d1, -0.2706056d1/), &
400 & C2 = (/0.0d0, 0.42981d-1, -0.293762d0, -0.1752461d1, 0.5682633d1, -0.3582633d1/)
402 REAL(8), DIMENSION(4), PARAMETER :: &
403 & C3 = (/0.5440d0, -0.39978d0, 0.25054d-1, -0.6714d-3/), &
404 & C4 = (/0.13822d1, -0.77857d0, 0.62767d-1, -0.20322d-2/), &
405 & C5 = (/-0.15861d1, -0.31082d0, -0.83751d-1, 0.38915d-2/)
407 REAL(8), DIMENSION(3), PARAMETER :: C6 = (/-0.4803d0, -0.82676d-1, 0.30302d-2/)
408 REAL(8), DIMENSION(2), PARAMETER :: C7 = (/0.164d0, 0.533d0/), &
409 & C8 = (/0.1736d0, 0.315d0/), &
410 & C9 = (/0.256d0, -0.635d-2/)
411 REAL(8), DIMENSION(2), PARAMETER :: G = (/-0.2273d1, 0.459d0/)
413 INTEGER(4) :: N2, i, i1
414 REAL(8) :: W, W1, pW, gamma, var1, var2, u, rsn, fac
415 REAL(8) :: Ndble, N25, summ2, ssumm2, a1, a2, range, sx
416 REAL(8) :: ssa, ssx, sax, ssassx, pi, mean, sd
417 REAL(8), DIMENSION(Nobs) :: X, XX, XSX, ASA
418 REAL(8), ALLOCATABLE :: a(:)
420 ! cxml function
421 REAL(8) :: dnrsq ! sum of squares
423 ! fortran(3f) routine
424 REAL(8), EXTERNAL :: derfc
425 ! ------------------------------------------------------------------------------
426 ! Shapiro-Wilk test (W test) for normality
427 ! ------------------------------------------------------------------------------
428 ! Based on the code AS R94 (P. Royston, Appl. Statist. 1995, 44, 547-551), obtained
429 ! from StatLib (http://lib.stat.cmu.edu/).
430 ! We implement only the case of an uncensored sample.
431 Tpass = .FALSE.
432 WRITE(*,*) '-----------------------------------------------------------------'
433 WRITE(*,*) '----------- Shapiro-Wilk test (W test) for normality -----------'
434 WRITE(*,*) '-----------------------------------------------------------------'
435 WRITE(*,*)
436 IF (Nobs < 3) THEN
437 WRITE(*,*) 'ERROR: The number of observations is < 3.'
438 WRITE(*,*) ' Unable to perform test.'
439 WRITE(*,*)
440 RETURN
441 END IF
442 IF (Nobs > 5000) THEN
443 WRITE(*,*) 'WARNING: The number of observations is > 5000.'
444 WRITE(*,*) ' The approximation for the p value may be inaccurate!'
445 WRITE(*,*)
446 END IF
448 ! Sort data in ascending order (cxml routine)
449 X = sample
450 CALL dsortq('A', Nobs, X, 1)
452 ! Check for zero range
453 range = X(Nobs) - X(1)
454 IF (range < small) THEN
455 WRITE(*,*) 'ERROR: Range of data is too small.'
456 WRITE(*,*)
457 RETURN
458 END IF
460 WRITE(*,*) 'Null hypothesis H0: The sample is drawn from a normal distribution.'
461 WRITE(*,*)
462 WRITE(*,*) 'Alternate hypothesis H1: The population is not normal.'
463 WRITE(*,*)
464 WRITE(*,*) 'Assumption: The sample is uncensored and untruncated.'
465 WRITE(*,*)
467 ! Calculate weights
468 Ndble = DBLE(Nobs)
469 N2 = Nobs/2 ! Integer division
470 ALLOCATE(a(N2))
472 IF (Nobs == 3) THEN
473 a(1) = 1.d0/SQRT(2.d0)
474 ELSE
475 N25 = Ndble + 0.25d0
476 DO i = 1, N2
477 a(i) = ppnd((DBLE(i) - 0.375d0)/N25)
478 END DO
480 summ2 = 2.d0*dnrsq(N2, a, 1)
481 ssumm2 = SQRT(summ2)
482 rsn = 1.d0 / SQRT(Ndble)
483 a1 = poly(C1, 6, rsn) - a(1)/ssumm2
485 ! Normalize coefficients
486 IF (Nobs > 5) THEN
487 i1 = 3
488 a2 = -a(2)/ssumm2 + poly(C2, 6, rsn)
489 fac = SQRT((summ2 - 2.d0*a(1)**2 - 2.d0*a(2)**2) &
490 & /(1.d0 - 2.d0*a1**2 - 2.d0*a2**2))
492 a(1) = a1
493 a(2) = a2
494 ELSE !Nobs = 4 or 5
495 i1 = 2
496 fac = sqrt((summ2 - 2.d0*a(1)**2)/(1.d0 - 2.d0*a1**2))
497 a(1) = a1
498 END IF
500 a(i1:) = -a(i1:)/fac
501 END IF
503 ! Scale with range
504 XX = X/range
506 ! Calculate W statistic as squared correlation between data and coefficients
507 ASA(1:N2) = -a
508 IF (2*N2.NE.Nobs) THEN
509 ASA(N2+1) = 0.d0
510 ASA(N2+2:Nobs) = a(N2:1:-1)
511 ELSE
512 ASA(N2+1:Nobs) = a(N2:1:-1)
513 END IF
515 sx = SUM(XX)/Ndble
516 XSX = XX - sx
517 ssx = dnrsq(Nobs, XSX, 1)
518 ssa = dnrsq(Nobs, ASA, 1)
519 sax = DOT_PRODUCT(XSX, ASA)
521 ! W1 equals (1-W) calculated to avoid excessive rounding error
522 ! for W very near 1 (a potential problem in very large samples)
523 ssassx = SQRT(ssa*ssx)
524 W1 = (ssassx - sax)*(ssassx + sax)/(ssa*ssx)
525 W = 1.d0 - W1
527 WRITE(*,'((A),F8.4)') ' Test statistic W: ', W
528 WRITE(*,*) '(W approaches 1 for a normal distribution)'
529 WRITE(*,*)
531 ! Calculate significance level for W (exact for N=3)
532 IF (Nobs == 3) THEN
533 pi = 4.d0*ATAN(1.d0)
534 pW = 6.d0/pi * (ASIN(SQRT(W)) - pi/3.d0)
535 ELSE
536 var1 = LOG(W1)
537 IF (Nobs <= 11) THEN
538 gamma = poly(G, 2, Ndble)
539 IF (var1 >= gamma) THEN
540 pW = small
541 ELSE
542 var1 = -LOG(gamma - var1)
543 mean = poly(C3, 4, Ndble)
544 sd = EXP(poly(C4, 4, Ndble))
545 u = (var1 - mean)/sd
546 pW = derfc(u/SQRT(2.d0))/2.d0
547 END IF
548 ELSE
549 var2 = LOG(Ndble)
550 mean = poly(C5, 4, var2)
551 sd = EXP(poly(C6, 3, var2))
552 u = (var1 - mean)/sd
553 pW = derfc(u/SQRT(2.d0))/2.d0
554 END IF
555 END IF
557 WRITE(*,'((A),E12.4)') ' Significance level of W: ', pW
558 WRITE(*,*)
559 WRITE(*,*) 'At significance level alpha = 0.05 (one-tailed test)'
560 IF (pW < 5.d-2) THEN
561 WRITE(*,*) '*** H0 IS REJECTED (no normal distribution) ***'
562 WRITE(*,*) 'Recommendation: Increase bin width.'
563 ELSE
564 WRITE(*,*) '*** H0 IS NOT REJECTED (normal distribution) ***'
565 Tpass = .TRUE.
566 END IF
567 WRITE(*,*) '-----------------------------------------------------------------'
569 DEALLOCATE(a)
570 RETURN
572 CONTAINS
573 ! ------------------------------------------------------------------------------
574 REAL(8) FUNCTION ppnd(p)
575 ! ------------------------------------------------------------------------------
576 ! Produces the normal deviate corresponding to given lower-tail area p, i.e.,
577 ! returns the value of z for which the cumulative density function of the
578 ! standard normal distribution is p. It is the inverse of the standard normal CDF.
579 ! z is accurate to about 1 part in 1.E16!
580 !
581 ! p = Phi(z) = 1/Sqrt[2*Pi]*Integrate[Exp[-t^2/2], {t, -Infinity, z}]
582 !
583 ! Algorithm AS 241, M.J. Wichura, Appl. Statist. 1988, 37, 477.
584 ! ------------------------------------------------------------------------------
585 IMPLICIT NONE
587 REAL(8), INTENT(in) :: p
589 REAL(8), PARAMETER :: ZERO = 0.D0, ONE = 1.D0, HALF = 0.5D0
590 REAL(8), PARAMETER :: SPLIT1 = 0.425D0, SPLIT2 = 5.D0
591 REAL(8), PARAMETER :: CONST1 = 0.180625D0, CONST2 = 1.6D0
593 ! Coefficients for p close to 0.5
594 REAL(8), PARAMETER :: A0 = 3.3871328727963666080D0, &
595 & A1 = 1.3314166789178437745D+2, &
596 & A2 = 1.9715909503065514427D+3, &
597 & A3 = 1.3731693765509461125D+4, &
598 & A4 = 4.5921953931549871457D+4, &
599 & A5 = 6.7265770927008700853D+4, &
600 & A6 = 3.3430575583588128105D+4, &
601 & A7 = 2.5090809287301226727D+3, &
602 & B1 = 4.2313330701600911252D+1, &
603 & B2 = 6.8718700749205790830D+2, &
604 & B3 = 5.3941960214247511077D+3, &
605 & B4 = 2.1213794301586595867D+4, &
606 & B5 = 3.9307895800092710610D+4, &
607 & B6 = 2.8729085735721942674D+4, &
608 & B7 = 5.2264952788528545610D+3
610 ! Coefficients for p not close to 0, 0.5 or 1
611 REAL(8), PARAMETER :: C0 = 1.42343711074968357734D0, &
612 & C1 = 4.63033784615654529590D0, &
613 & C2 = 5.76949722146069140550D0, &
614 & C3 = 3.64784832476320460504D0, &
615 & C4 = 1.27045825245236838258D0, &
616 & C5 = 2.41780725177450611770D-1, &
617 & C6 = 2.27238449892691845833D-2, &
618 & C7 = 7.74545014278341407640D-4, &
619 & D1 = 2.05319162663775882187D0, &
620 & D2 = 1.67638483018380384940D0, &
621 & D3 = 6.89767334985100004550D-1, &
622 & D4 = 1.48103976427480074590D-1, &
623 & D5 = 1.51986665636164571966D-2, &
624 & D6 = 5.47593808499534494600D-4, &
625 & D7 = 1.05075007164441684324D-9
627 ! Coefficients for p near 0 or 1
628 REAL(8), PARAMETER :: E0 = 6.65790464350110377720D0, &
629 & E1 = 5.46378491116411436990D0, &
630 & E2 = 1.78482653991729133580D0, &
631 & E3 = 2.96560571828504891230D-1, &
632 & E4 = 2.65321895265761230930D-2, &
633 & E5 = 1.24266094738807843860D-3, &
634 & E6 = 2.71155556874348757815D-5, &
635 & E7 = 2.01033439929228813265D-7, &
636 & F1 = 5.99832206555887937690D-1, &
637 & F2 = 1.36929880922735805310D-1, &
638 & F3 = 1.48753612908506148525D-2, &
639 & F4 = 7.86869131145613259100D-4, &
640 & F5 = 1.84631831751005468180D-5, &
641 & F6 = 1.42151175831644588870D-7, &
642 & F7 = 2.04426310338993978564D-15
644 REAL(8) :: q, r
645 ! ------------------------------------------------------------------------------
646 Q = P - HALF
647 IF (ABS(Q).LE.SPLIT1) THEN
648 R = CONST1 - Q**2
649 PPND = Q*(((((((A7*R + A6)*R + A5)*R + A4)*R + A3)*R + A2)*R + A1)*R + A0) &
650 & /(((((((B7*R + B6)*R + B5)*R + B4)*R + B3)*R + B2)*R + B1)*R + ONE)
651 RETURN
652 ELSE
653 IF (Q.LT.ZERO) THEN
654 R = P
655 ELSE
656 R = ONE - P
657 END IF
658 IF (R.LE.ZERO) THEN
659 STOP 'ERROR in function ppnd: p not in [0, 1].'
660 RETURN
661 END IF
662 R = SQRT(-LOG(R))
663 IF (R.LE.SPLIT2) THEN
664 R = R - CONST2
665 PPND = (((((((C7*R + C6)*R + C5)*R + C4)*R + C3)*R + C2)*R + C1)*R + C0) &
666 & /(((((((D7*R + D6)*R + D5)*R + D4)*R + D3)*R + D2)*R + D1)*R + ONE)
667 ELSE
668 R = R - SPLIT2
669 PPND = (((((((E7*R + E6)*R + E5)*R + E4)*R + E3)*R + E2)*R + E1)*R + E0) &
670 & /(((((((F7*R + F6)*R + F5)*R + F4)*R + F3)*R + F2)*R + F1)*R + ONE)
671 END IF
672 IF (Q.LT.ZERO) PPND = -PPND
673 RETURN
674 END IF
676 END FUNCTION PPND
678 ! ------------------------------------------------------------------------------
679 REAL(8) FUNCTION poly(c, nord, x)
680 ! ------------------------------------------------------------------------------
681 ! Calculates the algebraic polynomial of order nord-1 with
682 ! array of coefficients c. Zero-order coefficient is c(1)
683 !
684 ! Algorithm AS 181.2, J.P. Royston, Appl. Statist. 1982, 31, 176.
685 ! ------------------------------------------------------------------------------
686 IMPLICIT NONE
688 INTEGER(4), INTENT(in) :: nord
689 REAL(8), INTENT(in) :: c(nord)
690 REAL(8), INTENT(in) :: x
692 INTEGER(4) :: j
693 REAL(8) :: p
694 ! ------------------------------------------------------------------------------
695 IF (nord <= 0) STOP 'ERROR in function poly: nord <= 0.'
697 poly = c(1)
698 IF (nord == 1) RETURN
700 p = x*c(nord)
702 IF (nord == 2) THEN
703 poly = poly + p
704 RETURN
705 END IF
707 DO j = nord - 1, 2, -1
708 p = (p + c(j))*x
709 END DO
711 poly = poly + p
713 RETURN
714 END FUNCTION poly
716 END SUBROUTINE Shapiro_Wilk
719 !==============================================================================
720 SUBROUTINE Neumann(sample, var, Nobs, Tpass)
721 !==============================================================================
722 IMPLICIT NONE
724 INTEGER(4), INTENT(in) :: Nobs
725 REAL(8), INTENT(in) :: sample(Nobs), var
727 LOGICAL(4), INTENT(out) :: Tpass
729 REAL(8) :: sq, r, SDr, Ndble, u
730 REAL(8) :: diff(Nobs - 1)
732 ! cxml function
733 REAL(8) :: dnrsq ! sum of squares
736 ! Quantiles of standard normal distribution
737 REAL(8), PARAMETER :: u90 = 1.2816d0, &
738 & u95 = 1.6449d0, &
739 & u975 = 1.96d0, &
740 & u99 = 2.3263d0, &
741 & u995 = 2.5758d0
742 ! ------------------------------------------------------------------------------
743 ! von Neumann test for serial correlation
744 ! ------------------------------------------------------------------------------
745 Tpass = .FALSE.
746 WRITE(*,*) '-----------------------------------------------------------------'
747 WRITE(*,*) '---------- von Neumann test for serial correlation ----------'
748 WRITE(*,*) '-----------------------------------------------------------------'
749 WRITE(*,*)
750 IF (Nobs < 20) THEN
751 WRITE(*,*) 'WARNING: The number of observations is <20.'
752 WRITE(*,*) ' The approximation of normal distribution for r is inaccurate!'
753 WRITE(*,*)
754 END IF
755 WRITE(*,*) 'Null hypothesis H0: There is no positive correlation between'
756 WRITE(*,*) ' adjacent observations.'
757 WRITE(*,*)
758 WRITE(*,*) 'Alternate hypothesis H1: There is positive correlation.'
759 WRITE(*,*)
760 WRITE(*,*) 'The test statistic r is approx. normally distributed with mean = 1'
761 WRITE(*,*) 'and standard deviation SD(r).'
762 WRITE(*,*) 'Assumption: Sample is drawn from a normal distribution.'
763 WRITE(*,*)
765 Ndble = DBLE(Nobs)
767 ! Calculate mean-square successive difference
768 diff = sample(2:Nobs) - sample(1:Nobs - 1)
769 sq = dnrsq(Nobs - 1, diff, 1) / (2.d0*(Ndble - 1.d0))
771 ! Calculate test statistic r and its SD
772 r = sq / var
773 SDr = SQRT(Ndble/(Ndble - 1.d0))
775 ! Standard normal deviate of r
776 u = (r - 1.d0) / SDr
778 WRITE(*,'((A),F10.4)') ' Test statistic r: ', r
779 WRITE(*,'((A),F10.4)') ' Standard deviation SD(r): ', SDr
780 WRITE(*,'((A),F10.4)') ' Standard normal deviate u = S/SD(r): ', u
781 WRITE(*,*) '(r lies u SD''s from the expected mean 1)'
782 WRITE(*,*)
784 WRITE(*,*) 'At significance level alpha = 0.05 (one-tailed test)'
785 IF (ABS(u) > u95) THEN
786 WRITE(*,*) '*** H0 IS REJECTED (there is positive correlation) ***'
787 WRITE(*,*) 'Recommendation: Increase bin width.'
788 ELSE
789 WRITE(*,*) '*** H0 IS NOT REJECTED (no correlation) ***'
790 Tpass = .TRUE.
791 END IF
792 WRITE(*,*)
793 WRITE(*,*) 'Quantiles of standard normal distribution:'
794 WRITE(*,*) 'alpha u(1-alpha) u(1-alpha/2)'
795 WRITE(*,'(F6.2,2F10.4)') 0.1, u90, u95
796 WRITE(*,'(F6.2,2F10.4)') 0.05, u95, u975
797 WRITE(*,'(F6.2,2F10.4)') 0.01, u99, u995
798 WRITE(*,*) '-----------------------------------------------------------------'
800 RETURN
801 END SUBROUTINE Neumann
