Line data Source code
1 : !> \file mo_percentile.f90
2 : !> \brief \copybrief mo_percentile
3 : !> \details \copydetails mo_percentile
4 :
5 : !> \brief Median and percentiles.
6 : !> \details This module provides routines for median and percentiles.
7 : !> \changelog
8 : !! - Matthias Cuntz, Mar 2011
9 : !! - written
10 : !! - Juliane Mai, Jul 2012
11 : !! - different interpolation schemes in percentiles
12 : !! - Matthias Cuntz, Juliane Mai, Jul 2012
13 : !! - uses previous of ksmallest to half execution time
14 : !! - Matthias Cuntz, May 2014
15 : !! - removed numerical recipes
16 : !> \author Mathias Cuntz
17 : !> \date Mar 2011
18 : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
19 : !! FORCES is released under the LGPLv3+ license \license_note
20 : MODULE mo_percentile
21 :
22 : USE mo_kind, ONLY : i4, sp, dp
23 :
24 : Implicit NONE
25 :
26 : PUBLIC :: median ! Median
27 : PUBLIC :: n_element ! The n-th smallest value in an array
28 : PUBLIC :: percentile ! The value below which a certain percent of the input fall
29 : PUBLIC :: qmedian ! Quick median calculation, rearranges input
30 :
31 : ! ------------------------------------------------------------------
32 :
33 : !> \brief Median.
34 :
35 : !> \details
36 : !! Returns the median of the values in an array.
37 : !! If size is even, then the mean of the size/2 and size/2+1 element is the median.
38 : !!
39 : !! If an optinal mask is given, values only on those locations that correspond
40 : !! to true values in the mask are used.
41 : !!
42 : !! \b Example
43 : !!
44 : !! \code{.f90}
45 : !! vec = (/ 1.,2.,3.,4.,5.,6.,7.,8.,9.,10. /)
46 : !! ! Returns 5.5
47 : !! out = median(vec)
48 : !! \endcode
49 : !!
50 : !! See also example in test directory.
51 :
52 : !> \param[in] "real(sp/dp) :: vec(:)" 1D-array with input numbers
53 : !> \param[in] "logical, optional :: mask(:)" 1D-array of logical values with size(vec).
54 : !! If present, only those locations in vec
55 : !! corresponding to the true values in mask are used.
56 : !> \retval "real(sp/dp) :: out" Median of values in input array
57 :
58 : !> \author Matthias Cuntz
59 : !> \date Mar 2011
60 : !> \author Juliane Mai
61 : !> \date Jul 2012
62 : !! - uses previous of ksmallest to half execution time
63 :
64 : INTERFACE median
65 : MODULE PROCEDURE median_sp, median_dp
66 : END INTERFACE median
67 :
68 : ! ------------------------------------------------------------------
69 :
70 : !> \brief Nth smallest value in array.
71 :
72 : !> \details
73 : !! Returns the n-th smallest value in an array.
74 : !!
75 : !! If an optinal mask is given, values only on those locations that correspond
76 : !! to true values in the mask are used.
77 : !!
78 : !! \b Example
79 : !!
80 : !! \code{.f90}
81 : !! vec = (/ 1.,2.,3.,4.,5.,6.,7.,8.,9.,10. /)
82 : !! ! Returns 4
83 : !! out = n_element(vec,4)
84 : !! \endcode
85 : !!
86 : !! See also example in test directory.
87 : !!
88 : !! \b Literature
89 : !!
90 : !! 1. Niklaus Wirth. _"Algorithms and Data Structures"_. Prentice-Hall, Inc., 1985. ISBN 0-13-022005-1.
91 : !!
92 : !> \param[in] "real(sp/dp) :: vec(:)" 1D-array with input numbers
93 : !> \param[in] "integer(i4), optional :: n" Index of sorted array
94 : !> \param[in] "logical, optional :: mask(:)" 1D-array of logical values with size(vec).
95 : !! If present, only those locations in vec
96 : !! corresponding to the true values in mask are used.
97 : !> \param[out] "real(sp/dp) :: before" (n-1)-th smallest value in input array, e.g.
98 : !! for median/percentile calculations
99 : !> \param[out] "real(sp/dp) :: previous" Same as before
100 : !> \param[out] "real(sp/dp) :: after" (n+1)-th smallest value in input array
101 : !> \param[out] "real(sp/dp) :: next" Same as after
102 : !> \retval "real(sp/dp) :: out" N-th smallest value in input array
103 :
104 : !> \author Matthias Cuntz
105 : !> \date May 2014
106 : !! - based on qmedian
107 :
108 :
109 : INTERFACE n_element
110 : MODULE PROCEDURE n_element_dp, n_element_sp
111 : END INTERFACE n_element
112 :
113 : ! ------------------------------------------------------------------
114 :
115 : !> \brief Percentile.
116 :
117 : !> \details
118 : !! Returns the value below which a certain percent of array values fall.
119 : !!
120 : !! If an optinal mask is given, values only on those locations that correspond
121 : !! to true values in the mask are used.
122 : !!
123 : !! Different definitions can be applied to interpolate the stepwise CDF of the given data.
124 : !! 1. Inverse empirical CDF (no interpolation, default MATHEMATICA)
125 : !! 2. Linear interpolation (California method)
126 : !! 3. Element numbered closest
127 : !! 4. Linear interpolation (hydrologist method)
128 : !! 5. Mean-based estimate (Weibull method, default IMSL)
129 : !! 6. Mode-based estimate
130 : !! 7. Median-based estimate
131 : !! 8. normal distribution estimate
132 : !!
133 : !! See: http://reference.wolfram.com/mathematica/tutorial/BasicStatistics.html
134 : !!
135 : !! \b Example
136 : !!
137 : !! \code{.f90}
138 : !! vec = (/ 1.,2.,3.,4.,5.,6.,7.,8.,9.,10. /)
139 : !! ! Returns 10.
140 : !! out = percentile(vec,95.)
141 : !! ! Returns (10.,8)
142 : !! out = percentile(vec,(/95.,80./))
143 : !! \endcode
144 : !!
145 : !! See also example in test directory
146 :
147 : !> \param[in] "real(sp/dp) :: vec(:)" 1D-array with input numbers
148 : !> \param[in] "real(sp/dp) :: k[(:)]" Percentage of percentile, can be 1 dimensional
149 : !> \param[in] "logical, optional :: mask(:)" 1D-array of logical values with size(vec).
150 : !! If present, only those locations in vec
151 : !! corresponding to the true values in mask are used.
152 : !> \param[in] "integer(i4), optional :: mode_in" Specifies the interpolation scheme applied.\n
153 : !! Default:
154 : !! Inverse empirical CDF (no interpolation, default Mathematica)
155 : !! mode_in = 1_i4
156 : !> \retval "real(sp/dp) :: out[(size(k))]" k-th percentile of values in input array, can be
157 : !! 1 dimensional corresponding to k
158 :
159 : !> \author Matthias Cuntz
160 : !> \date Mar 2011
161 :
162 : !> \author Stephan Thober
163 : !> \date Dec 2011
164 : !! - added 1 dimensional version
165 :
166 : !> \author Juliane Mai
167 : !> \date Jul 2012
168 : !! - different interpolation schemes
169 :
170 : !> \authors Matthias Cuntz, Juliane Mai
171 : !> \date Jul 2012
172 : !! - uses previous of ksmallest to half execution time
173 :
174 : INTERFACE percentile
175 : MODULE PROCEDURE percentile_0d_sp, percentile_0d_dp, percentile_1d_sp, percentile_1d_dp
176 : END INTERFACE percentile
177 :
178 : ! ------------------------------------------------------------------
179 :
180 : !> \brief Quick median calculation
181 :
182 : !> \details
183 : !! Quick calculation of the median thereby rearranging the input array.
184 : !!
185 : !! \b Example
186 : !!
187 : !! \code{.f90}
188 : !! vec = (/ 1.,2.,3.,4.,5.,6.,7.,8.,9.,10. /)
189 : !! ! Returns 5.5
190 : !! out = qmedian(vec)
191 : !! \endcode
192 : !!
193 : !! See also example in test directory.
194 : !!
195 : !! \b Literature
196 : !! 1. Niklaus Wirth. _"Algorithms and Data Structures"_. Prentice-Hall, Inc., 1985. ISBN 0-13-022005-1.
197 : !!
198 : !> \param[inout] "real(sp/dp) :: vec(:)" 1D-array with input numbers.
199 : !! Will be rearranged on output.
200 : !> \retval "real(sp/dp) :: out" Median of values in input array
201 :
202 : !> \author Filip Hroch
203 : !! - as part of Munipack: http://munipack.physics.muni.cz
204 : !> \author Matthias Cuntz
205 : !> \date Jul 2012
206 : !! - function, k=n/2+1
207 : !! - real median for even n
208 :
209 : INTERFACE qmedian
210 : MODULE PROCEDURE qmedian_sp, qmedian_dp
211 : END INTERFACE qmedian
212 :
213 : ! ------------------------------------------------------------------
214 :
215 : PRIVATE
216 :
217 : ! ------------------------------------------------------------------
218 :
219 : CONTAINS
220 :
221 : ! ------------------------------------------------------------------
222 :
223 8 : FUNCTION median_dp(arrin, mask)
224 :
225 : IMPLICIT NONE
226 :
227 : REAL(dp), DIMENSION(:), INTENT(IN) :: arrin
228 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
229 : REAL(dp) :: median_dp
230 :
231 : INTEGER(i4) :: n
232 4 : REAL(dp), DIMENSION(:), ALLOCATABLE :: arr
233 4 : REAL(dp) :: tmp
234 :
235 4 : if (present(mask)) then
236 22 : n = count(mask)
237 6 : allocate(arr(n))
238 2 : arr = pack(arrin, mask)
239 :
240 2 : if (n < 2) stop 'median_dp: n < 2'
241 :
242 2 : if (mod(n, 2) == 0) then ! Even
243 1 : median_dp = n_element(arr, n / 2 + 1, previous = tmp)
244 1 : median_dp = 0.5_dp * (median_dp + tmp)
245 : else ! Odd
246 1 : median_dp = n_element(arr, (n + 1) / 2)
247 : end if
248 :
249 2 : deallocate(arr)
250 : else
251 2 : n = size(arrin)
252 2 : if (n < 2) stop 'median_dp: n < 2'
253 :
254 2 : if (mod(n, 2) == 0) then ! Even
255 1 : median_dp = n_element(arrin, n / 2 + 1, previous = tmp)
256 1 : median_dp = 0.5_dp * (median_dp + tmp)
257 : else ! Odd
258 1 : median_dp = n_element(arrin, (n + 1) / 2)
259 : end if
260 : end if
261 :
262 4 : END FUNCTION median_dp
263 :
264 :
265 8 : FUNCTION median_sp(arrin, mask)
266 :
267 : IMPLICIT NONE
268 :
269 : REAL(sp), DIMENSION(:), INTENT(IN) :: arrin
270 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
271 : REAL(sp) :: median_sp
272 :
273 : INTEGER(i4) :: n
274 4 : REAL(sp), DIMENSION(:), ALLOCATABLE :: arr
275 4 : REAL(sp) :: tmp
276 :
277 4 : if (present(mask)) then
278 22 : n = count(mask)
279 6 : allocate(arr(n))
280 2 : arr = pack(arrin, mask)
281 :
282 2 : if (n < 2) stop 'median_sp: n < 2'
283 :
284 2 : if (mod(n, 2) == 0) then ! Even
285 1 : median_sp = n_element(arr, n / 2 + 1, previous = tmp)
286 1 : median_sp = 0.5_sp * (median_sp + tmp)
287 : else ! Odd
288 1 : median_sp = n_element(arr, (n + 1) / 2)
289 : end if
290 :
291 2 : deallocate(arr)
292 : else
293 2 : n = size(arrin)
294 2 : if (n < 2) stop 'median_sp: n < 2'
295 :
296 2 : if (mod(n, 2) == 0) then ! Even
297 1 : median_sp = n_element(arrin, n / 2 + 1, previous = tmp)
298 1 : median_sp = 0.5_sp * (median_sp + tmp)
299 : else ! Odd
300 1 : median_sp = n_element(arrin, (n + 1) / 2)
301 : end if
302 : end if
303 :
304 8 : END FUNCTION median_sp
305 :
306 : ! ------------------------------------------------------------------
307 :
308 1580 : function n_element_dp(idat, n, mask, before, after, previous, next)
309 :
310 : IMPLICIT NONE
311 :
312 : real(dp), dimension(:), intent(in) :: idat
313 : integer(i4), intent(in) :: n
314 : logical, dimension(:), optional, intent(in) :: mask
315 : real(dp), optional, intent(out) :: before
316 : real(dp), optional, intent(out) :: after
317 : real(dp), optional, intent(out) :: previous
318 : real(dp), optional, intent(out) :: next
319 : real(dp) :: n_element_dp
320 :
321 790 : real(dp), dimension(:), allocatable :: dat
322 790 : real(dp) :: w
323 : integer(i4) :: nn, k
324 : integer(i4) :: l, r, i, j
325 :
326 790 : if (present(mask)) then
327 11 : nn = count(mask)
328 3 : allocate(dat(nn))
329 1 : dat = pack(idat, mask)
330 : else
331 789 : nn = size(idat)
332 2367 : allocate(dat(nn))
333 4531583 : dat = idat
334 : end if
335 :
336 790 : if (n < 1) stop 'n_element_dp: n < 1'
337 790 : if (n > nn) stop 'n_element_dp: n > size(pack(dat,mask))'
338 :
339 : !dat = idat
340 790 : nn = size(dat)
341 790 : k = n !nn/2 + 1
342 790 : l = 1
343 790 : r = nn
344 10775 : do while(l < r)
345 9985 : n_element_dp = dat(k)
346 9985 : i = l
347 9985 : j = r
348 : do
349 8735803 : do while(dat(i) < n_element_dp)
350 8735803 : i = i + 1
351 : enddo
352 4659710 : do while(n_element_dp < dat(j))
353 2271176 : j = j - 1
354 : enddo
355 2388534 : if (i <= j) then
356 2380140 : w = dat(i)
357 2380140 : dat(i) = dat(j)
358 2380140 : dat(j) = w
359 2380140 : i = i + 1
360 2380140 : j = j - 1
361 : end if
362 2388534 : if (i > j) exit
363 : enddo
364 9985 : if (j < k) l = i
365 9985 : if (k < i) r = j
366 : enddo
367 : ! if (mod(nn,2) == 0) then
368 : ! n_element_dp = 0.5_dp*(dat(k) + maxval(dat(:k-1)))
369 : ! else
370 : ! n_element_dp = dat(k)
371 : ! end if
372 790 : n_element_dp = dat(k)
373 :
374 790 : if (present(before)) before = maxval(dat(: k - 1))
375 908379 : if (present(previous)) previous = maxval(dat(: k - 1))
376 790 : if (present(after)) after = minval(dat(k + 1 :))
377 790 : if (present(next)) next = minval(dat(k + 1 :))
378 :
379 790 : deallocate(dat)
380 :
381 4 : end function n_element_dp
382 :
383 68 : function n_element_sp(idat, n, mask, before, after, previous, next)
384 :
385 : IMPLICIT NONE
386 :
387 : real(sp), dimension(:), intent(in) :: idat
388 : integer(i4), intent(in) :: n
389 : logical, dimension(:), optional, intent(in) :: mask
390 : real(sp), optional, intent(out) :: before
391 : real(sp), optional, intent(out) :: after
392 : real(sp), optional, intent(out) :: previous
393 : real(sp), optional, intent(out) :: next
394 : real(sp) :: n_element_sp
395 :
396 34 : real(sp), dimension(:), allocatable :: dat
397 34 : real(sp) :: w
398 : integer(i4) :: nn, k
399 : integer(i4) :: l, r, i, j
400 :
401 34 : if (present(mask)) then
402 11 : nn = count(mask)
403 3 : allocate(dat(nn))
404 1 : dat = pack(idat, mask)
405 : else
406 33 : nn = size(idat)
407 99 : allocate(dat(nn))
408 391 : dat = idat
409 : end if
410 :
411 34 : if (n < 1) stop 'n_element_sp: n < 1'
412 34 : if (n > nn) stop 'n_element_sp: n > size(pack(dat,mask))'
413 :
414 : !dat = idat
415 34 : nn = size(dat)
416 34 : k = n !nn/2 + 1
417 34 : l = 1
418 34 : r = nn
419 68 : do while(l < r)
420 34 : n_element_sp = dat(k)
421 34 : i = l
422 34 : j = r
423 : do
424 237 : do while(dat(i) < n_element_sp)
425 237 : i = i + 1
426 : enddo
427 131 : do while(n_element_sp < dat(j))
428 97 : j = j - 1
429 : enddo
430 34 : if (i <= j) then
431 34 : w = dat(i)
432 34 : dat(i) = dat(j)
433 34 : dat(j) = w
434 34 : i = i + 1
435 34 : j = j - 1
436 : end if
437 34 : if (i > j) exit
438 : enddo
439 34 : if (j < k) l = i
440 34 : if (k < i) r = j
441 : enddo
442 : ! if (mod(nn,2) == 0) then
443 : ! n_element_sp = 0.5_sp*(dat(k) + maxval(dat(:k-1)))
444 : ! else
445 : ! n_element_sp = dat(k)
446 : ! end if
447 34 : n_element_sp = dat(k)
448 :
449 34 : if (present(before)) before = maxval(dat(: k - 1))
450 178 : if (present(previous)) previous = maxval(dat(: k - 1))
451 34 : if (present(after)) after = minval(dat(k + 1 :))
452 34 : if (present(next)) next = minval(dat(k + 1 :))
453 :
454 34 : deallocate(dat)
455 :
456 790 : end function n_element_sp
457 :
458 : ! ------------------------------------------------------------------
459 :
460 1532 : FUNCTION percentile_0d_dp(arrin, k, mask, mode_in)
461 :
462 : IMPLICIT NONE
463 :
464 : REAL(dp), DIMENSION(:), INTENT(IN) :: arrin
465 : REAL(dp), INTENT(IN) :: k
466 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
467 : INTEGER(i4), OPTIONAL, INTENT(IN) :: mode_in
468 : REAL(dp) :: percentile_0d_dp
469 :
470 : INTEGER(i4) :: n, nn1, nn2
471 : INTEGER(i4) :: mode
472 766 : REAL(dp) :: kk, ks1, ks2
473 766 : REAL(dp), DIMENSION(:), ALLOCATABLE :: arr
474 :
475 766 : if (present(mask)) then
476 5497664 : n = count(mask)
477 : else
478 12 : n = size(arrin)
479 : end if
480 :
481 766 : if (present(mode_in)) then
482 764 : mode = mode_in
483 : else
484 : ! Default : Inverse empirical CDF
485 : mode = 1_i4
486 : end if
487 :
488 766 : if (n < 2) stop 'percentile_0d_dp: n < 2'
489 :
490 3 : select case (mode)
491 : ! Inverse empirical CDF: Mathematica default
492 : case(1_i4)
493 3 : kk = k / 100._dp * real(n, dp)
494 3 : nn1 = min(n, max(1_i4, ceiling(kk, kind = i4)))
495 3 : nn2 = nn1
496 :
497 : ! Linear interpolation (California method)
498 : case(2_i4)
499 1 : kk = k / 100._dp * real(n, dp)
500 1 : nn1 = min(n, max(1_i4, floor(kk, kind = i4)))
501 1 : nn2 = min(n, max(1_i4, ceiling(kk, kind = i4)))
502 :
503 : ! Element numbered closest
504 : case(3_i4)
505 1 : kk = 0.5_dp + k / 100._dp * real(n, dp)
506 1 : nn1 = min(n, max(1_i4, floor(kk, kind = i4)))
507 1 : nn2 = nn1
508 :
509 : ! Linear interpolation (hydrologist method)
510 : case(4_i4)
511 757 : kk = 0.5_dp + k / 100._dp * (real(n, dp))
512 757 : nn1 = min(n, max(1_i4, floor(kk, kind = i4)))
513 757 : nn2 = min(n, max(1_i4, ceiling(kk, kind = i4)))
514 :
515 : ! Mean-based estimate (Weibull method): IMSL default
516 : case(5_i4)
517 1 : kk = k / 100._dp * (real(n, dp) + 1._dp)
518 1 : nn1 = min(n, max(1_i4, floor(kk, kind = i4)))
519 1 : nn2 = min(n, max(1_i4, ceiling(kk, kind = i4)))
520 :
521 : ! Mode-based estimate
522 : case(6_i4)
523 1 : kk = 1.0_dp + k / 100._dp * (real(n, dp) - 1._dp)
524 1 : nn1 = min(n, max(1_i4, floor(kk, kind = i4)))
525 1 : nn2 = min(n, max(1_i4, ceiling(kk, kind = i4)))
526 :
527 : ! Median-based estimate
528 : case(7_i4)
529 1 : kk = 1.0_dp / 3.0_dp + k / 100._dp * (real(n, dp) + 1.0_dp / 3.0_dp)
530 1 : nn1 = min(n, max(1_i4, floor(kk, kind = i4)))
531 1 : nn2 = min(n, max(1_i4, ceiling(kk, kind = i4)))
532 :
533 : ! Normal distribution estimate
534 : case(8_i4)
535 1 : kk = 3.0_dp / 8.0_dp + k / 100._dp * (real(n, dp) + 1.0_dp / 4.0_dp)
536 1 : nn1 = min(n, max(1_i4, floor(kk, kind = i4)))
537 1 : nn2 = min(n, max(1_i4, ceiling(kk, kind = i4)))
538 :
539 : ! No valid mode
540 : case default
541 766 : stop 'percentile_0d_dp: mode > 8 not implemented'
542 :
543 : end select
544 :
545 766 : if (present(mask)) then
546 2262 : allocate(arr(n))
547 754 : arr = pack(arrin, mask)
548 754 : if (nn1 .eq. nn2) then
549 : ! no interpolation
550 1 : percentile_0d_dp = n_element(arr, nn1)
551 : else
552 : ! interpolation
553 753 : ks2 = n_element(arr, nn2, previous = ks1)
554 753 : percentile_0d_dp = ks1 + (ks2 - ks1) * (kk - real(nn1, dp))
555 : end if
556 754 : deallocate(arr)
557 : else
558 12 : if (nn1 .eq. nn2) then
559 : ! no interpolation
560 4 : percentile_0d_dp = n_element(arrin, nn1)
561 : else
562 : ! interpolation
563 8 : ks2 = n_element(arrin, nn2, previous = ks1)
564 8 : percentile_0d_dp = ks1 + (ks2 - ks1) * (kk - real(nn1, dp))
565 : end if
566 : end if
567 :
568 800 : END FUNCTION percentile_0d_dp
569 :
570 :
571 20 : FUNCTION percentile_0d_sp(arrin, k, mask, mode_in)
572 :
573 : IMPLICIT NONE
574 :
575 : REAL(sp), DIMENSION(:), INTENT(IN) :: arrin
576 : REAL(sp), INTENT(IN) :: k
577 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
578 : INTEGER(i4), OPTIONAL, INTENT(IN) :: mode_in
579 : REAL(sp) :: percentile_0d_sp
580 :
581 : INTEGER(i4) :: n, nn1, nn2
582 : INTEGER(i4) :: mode
583 10 : REAL(sp) :: kk, ks1, ks2
584 10 : REAL(sp), DIMENSION(:), ALLOCATABLE :: arr
585 :
586 10 : if (present(mask)) then
587 11 : n = count(mask)
588 : else
589 9 : n = size(arrin)
590 : end if
591 :
592 10 : if (present(mode_in)) then
593 8 : mode = mode_in
594 : else
595 : ! Default : Inverse empirical CDF
596 : mode = 1_i4
597 : end if
598 :
599 10 : if (n < 2) stop 'percentile_0d_sp: n < 2'
600 :
601 3 : select case (mode)
602 : ! Inverse empirical CDF: Mathematica default
603 : case(1_i4)
604 3 : kk = k / 100._sp * real(n, sp)
605 3 : nn1 = min(n, max(1_i4, ceiling(kk, kind = i4)))
606 3 : nn2 = nn1
607 :
608 : ! Linear interpolation (California method)
609 : case(2_i4)
610 1 : kk = k / 100._sp * real(n, sp)
611 1 : nn1 = min(n, max(1_i4, floor(kk, kind = i4)))
612 1 : nn2 = min(n, max(1_i4, ceiling(kk, kind = i4)))
613 :
614 : ! Element numbered closest
615 : case(3_i4)
616 1 : kk = 0.5_sp + k / 100._sp * real(n, sp)
617 1 : nn1 = min(n, max(1_i4, floor(kk, kind = i4)))
618 1 : nn2 = nn1
619 :
620 : ! Linear interpolation (hydrologist method)
621 : case(4_i4)
622 1 : kk = 0.5_sp + k / 100._sp * (real(n, sp))
623 1 : nn1 = min(n, max(1_i4, floor(kk, kind = i4)))
624 1 : nn2 = min(n, max(1_i4, ceiling(kk, kind = i4)))
625 :
626 : ! Mean-based estimate (Weibull method): IMSL default
627 : case(5_i4)
628 1 : kk = k / 100._sp * (real(n, sp) + 1._sp)
629 1 : nn1 = min(n, max(1_i4, floor(kk, kind = i4)))
630 1 : nn2 = min(n, max(1_i4, ceiling(kk, kind = i4)))
631 :
632 : ! Mode-based estimate
633 : case(6_i4)
634 1 : kk = 1.0_sp + k / 100._sp * (real(n, sp) - 1._sp)
635 1 : nn1 = min(n, max(1_i4, floor(kk, kind = i4)))
636 1 : nn2 = min(n, max(1_i4, ceiling(kk, kind = i4)))
637 :
638 : ! Median-based estimate
639 : case(7_i4)
640 1 : kk = 1.0_sp / 3.0_sp + k / 100._sp * (real(n, sp) + 1.0_sp / 3.0_sp)
641 1 : nn1 = min(n, max(1_i4, floor(kk, kind = i4)))
642 1 : nn2 = min(n, max(1_i4, ceiling(kk, kind = i4)))
643 :
644 : ! Normal distribution estimate
645 : case(8_i4)
646 1 : kk = 3.0_sp / 8.0_sp + k / 100._sp * (real(n, sp) + 1.0_sp / 4.0_sp)
647 1 : nn1 = min(n, max(1_i4, floor(kk, kind = i4)))
648 1 : nn2 = min(n, max(1_i4, ceiling(kk, kind = i4)))
649 :
650 : ! No valid mode
651 : case default
652 10 : stop 'percentile_0d_sp: mode > 8 not implemented'
653 :
654 : end select
655 :
656 10 : if (present(mask)) then
657 3 : allocate(arr(n))
658 1 : arr = pack(arrin, mask)
659 1 : if (nn1 .eq. nn2) then
660 : ! no interpolation
661 1 : percentile_0d_sp = n_element(arr, nn1)
662 : else
663 : ! interpolation
664 0 : ks2 = n_element(arr, nn2, previous = ks1)
665 0 : percentile_0d_sp = ks1 + (ks2 - ks1) * (kk - real(nn1, sp))
666 : end if
667 1 : deallocate(arr)
668 : else
669 9 : if (nn1 .eq. nn2) then
670 : ! no interpolation
671 4 : percentile_0d_sp = n_element(arrin, nn1)
672 : else
673 : ! interpolation
674 5 : ks2 = n_element(arrin, nn2, previous = ks1)
675 5 : percentile_0d_sp = ks1 + (ks2 - ks1) * (kk - real(nn1, sp))
676 : end if
677 : end if
678 :
679 776 : END FUNCTION percentile_0d_sp
680 :
681 :
682 37 : function percentile_1d_dp(arrin, k, mask, mode_in)
683 :
684 : IMPLICIT NONE
685 :
686 : REAL(dp), DIMENSION(:), INTENT(IN) :: arrin
687 : REAL(dp), DIMENSION(:), INTENT(IN) :: k
688 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
689 : INTEGER(i4), OPTIONAL, INTENT(IN) :: mode_in
690 :
691 : REAL(dp), DIMENSION(size(k)) :: percentile_1d_dp
692 :
693 : INTEGER(i4) :: i, n
694 : INTEGER(i4) :: mode
695 9 : INTEGER(i4), DIMENSION(size(k)) :: nn1, nn2
696 45 : REAL(dp), DIMENSION(size(k)) :: kk
697 9 : REAL(dp) :: ks1, ks2
698 9 : REAL(dp), DIMENSION(:), ALLOCATABLE :: arr
699 :
700 9 : if (present(mask)) then
701 0 : n = count(mask)
702 : else
703 9 : n = size(arrin)
704 : end if
705 :
706 9 : if (present(mode_in)) then
707 8 : mode = mode_in
708 : else
709 : ! Default : Inverse empirical CDF
710 : mode = 1_i4
711 : end if
712 :
713 : ! check consistency
714 : !if (size(k) > size(arr)) stop 'percentile_1d_dp: more Quantiles than data: size(k) > size(arr)'
715 9 : if (n < 2) stop 'percentile_1d_dp: n < 2'
716 :
717 2 : select case (mode)
718 : ! Inverse empirical CDF: Mathematica default
719 : case(1_i4)
720 6 : kk(:) = k(:) / 100._dp * real(n, dp)
721 6 : nn1(:) = min(n, max(1_i4, ceiling(kk(:), kind = i4)))
722 6 : nn2 = nn1
723 :
724 : ! Linear interpolation (California method)
725 : case(2_i4)
726 3 : kk(:) = k(:) / 100._dp * real(n, dp)
727 3 : nn1(:) = min(n, max(1_i4, floor(kk(:), kind = i4)))
728 3 : nn2(:) = min(n, max(1_i4, ceiling(kk(:), kind = i4)))
729 :
730 : ! Element numbered closest
731 : case(3_i4)
732 3 : kk(:) = 0.5_dp + k(:) / 100._dp * real(n, dp)
733 3 : nn1(:) = min(n, max(1_i4, floor(kk(:), kind = i4)))
734 3 : nn2 = nn1
735 :
736 : ! Linear interpolation (hydrologist method)
737 : case(4_i4)
738 3 : kk(:) = 0.5_dp + k(:) / 100._dp * (real(n, dp))
739 3 : nn1(:) = min(n, max(1_i4, floor(kk(:), kind = i4)))
740 3 : nn2(:) = min(n, max(1_i4, ceiling(kk(:), kind = i4)))
741 :
742 : ! Mean-based estimate (Weibull method): IMSL default
743 : case(5_i4)
744 3 : kk(:) = k(:) / 100._dp * (real(n, dp) + 1._dp)
745 3 : nn1(:) = min(n, max(1_i4, floor(kk(:), kind = i4)))
746 3 : nn2(:) = min(n, max(1_i4, ceiling(kk(:), kind = i4)))
747 :
748 : ! Mode-based estimate
749 : case(6_i4)
750 3 : kk(:) = 1.0_dp + k(:) / 100._dp * (real(n, dp) - 1._dp)
751 3 : nn1(:) = min(n, max(1_i4, floor(kk(:), kind = i4)))
752 3 : nn2(:) = min(n, max(1_i4, ceiling(kk(:), kind = i4)))
753 :
754 : ! Median-based estimate
755 : case(7_i4)
756 3 : kk(:) = 1.0_dp / 3.0_dp + k(:) / 100._dp * (real(n, dp) + 1.0_dp / 3.0_dp)
757 3 : nn1(:) = min(n, max(1_i4, floor(kk(:), kind = i4)))
758 3 : nn2(:) = min(n, max(1_i4, ceiling(kk(:), kind = i4)))
759 :
760 : ! Normal distribution estimate
761 : case(8_i4)
762 3 : kk(:) = 3.0_dp / 8.0_dp + k(:) / 100._dp * (real(n, dp) + 1.0_dp / 4.0_dp)
763 3 : nn1(:) = min(n, max(1_i4, floor(kk(:), kind = i4)))
764 3 : nn2(:) = min(n, max(1_i4, ceiling(kk(:), kind = i4)))
765 :
766 : ! No valid mode
767 : case default
768 9 : stop 'percentile_1d_dp: mode > 8 not implemented'
769 :
770 : end select
771 :
772 9 : if (present(mask)) then
773 0 : allocate(arr(n))
774 0 : arr = pack(arrin, mask)
775 0 : do i = 1, size(k)
776 0 : if (nn1(i) .eq. nn2(i)) then
777 : ! no interpolation
778 0 : percentile_1d_dp(i) = n_element(arr, nn1(i))
779 : else
780 : ! interpolation
781 0 : ks2 = n_element(arr, nn2(i), previous = ks1)
782 0 : percentile_1d_dp(i) = ks1 + (ks2 - ks1) * (kk(i) - real(nn1(i), dp))
783 : end if
784 : end do
785 0 : deallocate(arr)
786 : else
787 27 : do i = 1, size(k)
788 27 : if (nn1(i) .eq. nn2(i)) then
789 : ! no interpolation
790 8 : percentile_1d_dp(i) = n_element(arrin, nn1(i))
791 : else
792 : ! interpolation
793 10 : ks2 = n_element(arrin, nn2(i), previous = ks1)
794 10 : percentile_1d_dp(i) = ks1 + (ks2 - ks1) * (kk(i) - real(nn1(i), dp))
795 : end if
796 : end do
797 : end if
798 :
799 9 : END function percentile_1d_dp
800 :
801 :
802 36 : function percentile_1d_sp(arrin, k, mask, mode_in)
803 :
804 : IMPLICIT NONE
805 :
806 : REAL(sp), DIMENSION(:), INTENT(IN) :: arrin
807 : REAL(sp), DIMENSION(:), INTENT(IN) :: k
808 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
809 : INTEGER(i4), OPTIONAL, INTENT(IN) :: mode_in
810 :
811 : REAL(sp), DIMENSION(size(k)) :: percentile_1d_sp
812 :
813 : INTEGER(i4) :: i, n
814 : INTEGER(i4) :: mode
815 9 : INTEGER(i4), DIMENSION(size(k)) :: nn1, nn2
816 45 : REAL(sp), DIMENSION(size(k)) :: kk
817 9 : REAL(sp) :: ks1, ks2
818 9 : REAL(sp), DIMENSION(:), ALLOCATABLE :: arr
819 :
820 9 : if (present(mask)) then
821 0 : n = count(mask)
822 : else
823 9 : n = size(arrin)
824 : end if
825 :
826 9 : if (present(mode_in)) then
827 8 : mode = mode_in
828 : else
829 : ! Default : Inverse empirical CDF
830 : mode = 1_i4
831 : end if
832 :
833 : ! check consistency
834 : !if (size(k) > size(arr)) stop 'percentile_1d_sp: more Quantiles than data: size(k) > size(arr)'
835 9 : if (n < 2) stop 'percentile_1d_sp: n < 2'
836 :
837 2 : select case (mode)
838 : ! Inverse empirical CDF: Mathematica default
839 : case(1_i4)
840 6 : kk(:) = k(:) / 100._sp * real(n, sp)
841 6 : nn1(:) = min(n, max(1_i4, ceiling(kk(:), kind = i4)))
842 6 : nn2 = nn1
843 :
844 : ! Linear interpolation (California method)
845 : case(2_i4)
846 3 : kk(:) = k(:) / 100._sp * real(n, sp)
847 3 : nn1(:) = min(n, max(1_i4, floor(kk(:), kind = i4)))
848 3 : nn2(:) = min(n, max(1_i4, ceiling(kk(:), kind = i4)))
849 :
850 : ! Element numbered closest
851 : case(3_i4)
852 3 : kk(:) = 0.5_sp + k(:) / 100._sp * real(n, sp)
853 3 : nn1(:) = min(n, max(1_i4, floor(kk(:), kind = i4)))
854 3 : nn2 = nn1
855 :
856 : ! Linear interpolation (hydrologist method)
857 : case(4_i4)
858 3 : kk(:) = 0.5_sp + k(:) / 100._sp * (real(n, sp))
859 3 : nn1(:) = min(n, max(1_i4, floor(kk(:), kind = i4)))
860 3 : nn2(:) = min(n, max(1_i4, ceiling(kk(:), kind = i4)))
861 :
862 : ! Mean-based estimate (Weibull method): IMSL default
863 : case(5_i4)
864 3 : kk(:) = k(:) / 100._sp * (real(n, sp) + 1._sp)
865 3 : nn1(:) = min(n, max(1_i4, floor(kk(:), kind = i4)))
866 3 : nn2(:) = min(n, max(1_i4, ceiling(kk(:), kind = i4)))
867 :
868 : ! Mode-based estimate
869 : case(6_i4)
870 3 : kk(:) = 1.0_sp + k(:) / 100._sp * (real(n, sp) - 1._sp)
871 3 : nn1(:) = min(n, max(1_i4, floor(kk(:), kind = i4)))
872 3 : nn2(:) = min(n, max(1_i4, ceiling(kk(:), kind = i4)))
873 :
874 : ! Median-based estimate
875 : case(7_i4)
876 3 : kk(:) = 1.0_sp / 3.0_sp + k(:) / 100._sp * (real(n, sp) + 1.0_sp / 3.0_sp)
877 3 : nn1(:) = min(n, max(1_i4, floor(kk(:), kind = i4)))
878 3 : nn2(:) = min(n, max(1_i4, ceiling(kk(:), kind = i4)))
879 :
880 : ! Normal distribution estimate
881 : case(8_i4)
882 3 : kk(:) = 3.0_sp / 8.0_sp + k(:) / 100._sp * (real(n, sp) + 1.0_sp / 4.0_sp)
883 3 : nn1(:) = min(n, max(1_i4, floor(kk(:), kind = i4)))
884 3 : nn2(:) = min(n, max(1_i4, ceiling(kk(:), kind = i4)))
885 :
886 : ! No valid mode
887 : case default
888 9 : stop 'percentile_1d_sp: mode > 8 not implemented'
889 :
890 : end select
891 :
892 9 : if (present(mask)) then
893 0 : allocate(arr(n))
894 0 : arr = pack(arrin, mask)
895 0 : do i = 1, size(k)
896 0 : if (nn1(i) .eq. nn2(i)) then
897 : ! no interpolation
898 0 : percentile_1d_sp(i) = n_element(arr, nn1(i))
899 : else
900 : ! interpolation
901 0 : ks2 = n_element(arr, nn2(i), previous = ks1)
902 0 : percentile_1d_sp(i) = ks1 + (ks2 - ks1) * (kk(i) - real(nn1(i), sp))
903 : end if
904 : end do
905 0 : deallocate(arr)
906 : else
907 27 : do i = 1, size(k)
908 27 : if (nn1(i) .eq. nn2(i)) then
909 : ! no interpolation
910 8 : percentile_1d_sp(i) = n_element(arrin, nn1(i))
911 : else
912 : ! interpolation
913 10 : ks2 = n_element(arrin, nn2(i), previous = ks1)
914 10 : percentile_1d_sp(i) = ks1 + (ks2 - ks1) * (kk(i) - real(nn1(i), sp))
915 : end if
916 : end do
917 : end if
918 :
919 9 : END function percentile_1d_sp
920 :
921 : ! ------------------------------------------------------------------
922 :
923 2 : function qmedian_dp(dat)
924 :
925 : IMPLICIT NONE
926 :
927 : real(dp), dimension(:), intent(inout) :: dat
928 : real(dp) :: qmedian_dp
929 :
930 2 : real(dp) :: w
931 : integer(i4) :: n, k
932 : integer(i4) :: l, r, i, j
933 :
934 2 : n = size(dat)
935 2 : k = n / 2 + 1
936 2 : l = 1
937 2 : r = n
938 4 : do while(l < r)
939 2 : qmedian_dp = dat(k)
940 2 : i = l
941 2 : j = r
942 : do
943 11 : do while(dat(i) < qmedian_dp)
944 11 : i = i + 1
945 : enddo
946 10 : do while(qmedian_dp < dat(j))
947 8 : j = j - 1
948 : enddo
949 2 : if (i <= j) then
950 2 : w = dat(i)
951 2 : dat(i) = dat(j)
952 2 : dat(j) = w
953 2 : i = i + 1
954 2 : j = j - 1
955 : end if
956 2 : if (i > j) exit
957 : enddo
958 2 : if (j < k) l = i
959 2 : if (k < i) r = j
960 : enddo
961 2 : if (mod(n, 2) == 0) then
962 8 : qmedian_dp = 0.5_dp * (dat(k) + maxval(dat(: k - 1)))
963 : else
964 1 : qmedian_dp = dat(k)
965 : end if
966 :
967 9 : end function qmedian_dp
968 :
969 :
970 2 : function qmedian_sp(dat)
971 :
972 : IMPLICIT NONE
973 :
974 : real(sp), dimension(:), intent(inout) :: dat
975 : real(sp) :: qmedian_sp
976 :
977 2 : real(sp) :: w
978 : integer(i4) :: n, k
979 : integer(i4) :: l, r, i, j
980 :
981 2 : n = size(dat)
982 2 : k = n / 2 + 1
983 2 : l = 1
984 2 : r = n
985 4 : do while(l < r)
986 2 : qmedian_sp = dat(k)
987 2 : i = l
988 2 : j = r
989 : do
990 11 : do while(dat(i) < qmedian_sp)
991 11 : i = i + 1
992 : enddo
993 10 : do while(qmedian_sp < dat(j))
994 8 : j = j - 1
995 : enddo
996 2 : if (i <= j) then
997 2 : w = dat(i)
998 2 : dat(i) = dat(j)
999 2 : dat(j) = w
1000 2 : i = i + 1
1001 2 : j = j - 1
1002 : end if
1003 2 : if (i > j) exit
1004 : enddo
1005 2 : if (j < k) l = i
1006 2 : if (k < i) r = j
1007 : enddo
1008 2 : if (mod(n, 2) == 0) then
1009 8 : qmedian_sp = 0.5_sp * (dat(k) + maxval(dat(: k - 1)))
1010 : else
1011 1 : qmedian_sp = dat(k)
1012 : end if
1013 :
1014 2 : end function qmedian_sp
1015 :
1016 : ! ------------------------------------------------------------------
1017 :
1018 : END MODULE mo_percentile
|
