Line data Source code
1 : !> \file mo_moment.f90
2 : !> \brief \copybrief mo_moment
3 : !> \details \copydetails mo_moment
4 :
5 : !> \brief Statistical moments.
6 : !> \details This module contains routines for the masked calculation of
7 : !! statistical properties such as moments and mixed moments of input vectors
8 : !! \note all except variance and standard deviation are population and not sample moments,
9 : !! i.e. they are normally divided by n and not (n-1)
10 : !> \par Literature
11 : !! Central moments and error variances
12 : !! LH Benedict & RD Gould, Towards better uncertainty estimates for turbulence statistics,
13 : !! Experiments in Fluids 22, 129-136, 1996
14 : !> \changelog
15 : !! - Matthias Cuntz, Nov 2011
16 : !! - written
17 : !! - Matthias Cuntz, Dec 2011
18 : !! - mod. correlation, covariance
19 : !! - M. Schroen, Sep 2015
20 : !! - average/mean for single value
21 : !! - S. Mueller, Dec 2019
22 : !! - remove citations for common sence
23 : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
24 : !! FORCES is released under the LGPLv3+ license \license_note
25 : MODULE mo_moment
26 :
27 : USE mo_kind, ONLY : i4, sp, dp
28 :
29 : IMPLICIT NONE
30 :
31 : PUBLIC :: absdev ! Mean absolute deviation from mean of an array
32 : PUBLIC :: average ! 1st moment of an array (same as mean)
33 : PUBLIC :: central_moment ! Arbitrary central moment of an array
34 : PUBLIC :: central_moment_var ! Central moment error variance
35 : PUBLIC :: correlation ! Correlation between two arrays
36 : PUBLIC :: covariance ! Covariance between two arrays
37 : PUBLIC :: kurtosis ! 4th moment of an array
38 : PUBLIC :: mean ! 1st moment of an array
39 : PUBLIC :: mixed_central_moment ! Arbitrary mixed central moment
40 : PUBLIC :: mixed_central_moment_var ! Arbitrary mixed central moment error variance
41 : PUBLIC :: moment ! 1st to 4th moments of an array
42 : PUBLIC :: skewness ! 3rd moment of an array
43 : PUBLIC :: stddev ! Sqrt of 2nd moment of an array
44 : PUBLIC :: variance ! 2nd moment of an array
45 :
46 : ! ------------------------------------------------------------------
47 :
48 : !> \brief Mean absolute deviation from mean.
49 :
50 : !> \details
51 : !! Calculates the mean absolute deviations from the mean
52 : !!
53 : !! \f[ ABSDEV = \sum_i\frac{|x_i-\bar x|}{N} \f]
54 : !!
55 : !! If an optinal mask is given, the calculations are over those locations that correspond to true values in the mask.
56 : !! \f$ x\f$ can be single or double precision. The result will have the same numerical precision.\n
57 : !!
58 : !! \b Example
59 : !!
60 : !! vec = (/ 1., 2, 3., -999., 5., 6. /)
61 : !! m = absdev(vec, mask=(vec >= 0.))
62 : !!
63 : !! See also example in pf_tests directory.
64 :
65 :
66 : !> \param[in] "real(sp/dp) :: dat(:)" 1D-array with input numbers.
67 : !> \param[in] "logical, optional :: mask(:)" 1D-array of logical values with `size(dat)`.
68 : !! If present, only those locations in `vec`
69 : !! corresponding to the true values in mask are used.
70 : !> \retval "real(sp/dp) :: absdev" Mean absolute deviations from average.
71 :
72 : !> \note Input values must be floating points.
73 :
74 : !> \author Matthias Cuntz
75 : !> \date Nov 2011
76 :
77 : ! ------------------------------------------------------------------
78 :
79 : INTERFACE absdev
80 : MODULE PROCEDURE absdev_sp, absdev_dp
81 : END INTERFACE absdev
82 :
83 : ! ------------------------------------------------------------------
84 :
85 : !> \brief Mean of vector.
86 :
87 : !> \details
88 : !! Calculates the average value of a vector, i.e. the first moment of a series of numbers:
89 : !!
90 : !! \f[ AVE = \sum_i \frac{x_i}{N} \f]
91 : !!
92 : !! If an optional mask is given, the average is only over those locations that correspond to true values in the mask.
93 : !! \f$ x\f$ can be single or double precision. The result will have the same numerical precision.
94 : !!
95 : !! \b Example
96 : !!
97 : !! vec = (/ 1., 2, 3., -999., 5., 6. /)
98 : !! m = average(vec, mask=(vec >= 0.))
99 : !! --> m = 3.4
100 : !!
101 : !! See also example in pf_tests directory.
102 :
103 : !> \param[in] "real(sp/dp) :: dat(:)" 1D-array with input numbers
104 :
105 : !> \retval "real(sp/dp) :: average" Average of all elements in dat
106 :
107 :
108 : !> \param[in] "logical, optional :: mask(:)" 1D-array of logical values with `size(dat)`.
109 : !! If present, only those locations in dat
110 : !! corresponding to the true values in mask are used.
111 :
112 : !> \note Input values must be floating points.
113 :
114 : !> \author Matthias Cuntz
115 : !> \date Nov 2011
116 :
117 : ! ------------------------------------------------------------------
118 :
119 : INTERFACE average
120 : MODULE PROCEDURE average_sp, average_dp
121 : END INTERFACE average
122 :
123 : ! ------------------------------------------------------------------
124 :
125 : !> \brief R-central moment
126 :
127 : !> \details
128 : !! Calculates the central moment of a vector, i.e. the r-central moment of a series of numbers:
129 : !!
130 : !! \f[ \mu_r = \sum_i\frac{(x_i-\bar x)^r}{N} \f]
131 : !!
132 : !! Note that the variance is the second central moment: `variance(x) = central_moment(x,2)`\n
133 : !!
134 : !! If an optinal mask is given, the average is only over those locations that correspond to true values in the mask.
135 : !! x can be single or double precision. The result will have the same numerical precision.
136 : !!
137 : !! \b Example
138 : !!
139 : !! vec = (/ 1., 2, 3., 4., 5., 6. /)
140 : !! m = central_moment(vec, 2, mask=(vec >= 0.))
141 : !! --> m = 2.91666
142 : !!
143 : !! See also example in pf_tests directory.
144 : !!
145 : !> \b Literature
146 : !! 1. LH Benedict & RD Gould, _Towards better uncertainty estimates for turbulence statistics_.
147 : !! Experiments in Fluids 22, 129-136, 1996
148 : !!
149 : !> \param[in] "real(sp/dp) :: dat(:)" 1D-array with input numbers.
150 : !> \param[in] "integer(i4) :: r" Order of the central moment, i.e. r=2 is variance.
151 : !> \param[in] "logical, optional :: mask(:)" 1D-array of logical values with size(dat).
152 : !! If present, only those locations in `dat`
153 : !! corresponding to the true values in mask are used.
154 : !> \retval "real(sp/dp) :: central_moment" R-th central moment of elements in `dat`.
155 :
156 : !> \note Input values must be floating points.
157 :
158 : !> \author Matthias Cuntz
159 : !> \date Nov 2011
160 :
161 : ! ------------------------------------------------------------------
162 : INTERFACE central_moment
163 : MODULE PROCEDURE central_moment_sp, central_moment_dp
164 : END INTERFACE central_moment
165 :
166 : ! ------------------------------------------------------------------
167 :
168 : !> \brief R-central moment variance
169 :
170 : !> \details
171 : !! Calculates the sampling variance of the central moment of a vector.
172 : !! `central_moment_var` is something like the "error variance" of the r-th order central moment sampling statistic.\n
173 : !!
174 : !! If an optinal mask is given, the average is only over those locations that correspond to true values in the mask.
175 : !! \f$x\f$ can be single or double precision. The result will have the same numerical precision.\n
176 : !!
177 : !! \b Example
178 : !!
179 : !! vec = (/ 1., 2, 3., -999., 5., 6. /)
180 : !! m = central_moment(vec, 2, mask=(vec >= 0.))
181 : !! em = central_moment_var(vec, 2, mask=(vec >= 0.))
182 : !!
183 : !! See also example in pf_tests directory.
184 : !!
185 : !! \b Literature
186 : !! 1. LH Benedict & RD Gould, _Towards better uncertainty estimates for turbulence statistics_,
187 : !! Experiments in Fluids 22, 129-136, 1996
188 : !!
189 : !> \param[in] "real(sp/dp) :: dat(:)" 1D-array with input numbers.
190 : !> \param[in] "integer(i4) :: r" Order of the central moment, i.e. r=2 is variance.
191 : !> \param[in] "logical, optional :: mask(:)" 1D-array of logical values with `size(dat)`.
192 : !! If present, only those locations in dat
193 : !! corresponding to the true values in mask are used.
194 : !> \retval "real(sp/dp) :: central_moment_var" Sampling variance of r-th central moment of elements in dat
195 :
196 : !> \note Input values must be floating points.
197 :
198 : !> \author Matthias Cuntz
199 : !> \date Nov 2011
200 : INTERFACE central_moment_var
201 : MODULE PROCEDURE central_moment_var_sp, central_moment_var_dp
202 : END INTERFACE central_moment_var
203 :
204 : ! ------------------------------------------------------------------
205 :
206 : !> \brief Correlation between two vectors.
207 :
208 : !> \details
209 : !! Calculates the correlation between two input vectors, i.e. the covariance divided by the standard deviations:\n
210 : !!
211 : !! \f[\langle x y\rangle = \sum_i\frac{(x_i-\bar x)(y_i-\bar y)}{N\sqrt{\mu_2(x)\mu_2(y)}}\f]
212 : !!
213 : !! If an optinal mask is given, the average is only over those locations that correspond to true values in the mask.
214 : !! \f$x\f$ can be single or double precision. The result will have the same numerical precision.
215 : !!
216 : !! \b Example
217 : !!
218 : !! vec1 = (/ 1., 2, 3., -999., 5., 6. /)
219 : !! vec2 = (/ 1.3, 2.7, 3.9, 5.1, 6., 8. /)
220 : !! m = correlation(vec1, vec2, mask=((vec1 >= 0.) .and. (vec2 >= 0.)))
221 : !!
222 : !! See also example in pf_tests directory.
223 :
224 : !> \param[in] "real(sp/dp) :: x(:)" First 1D-array with input numbers.
225 : !> \param[in] "real(sp/dp) :: y(:)" Second 1D-array with input numbers.
226 : !> \retval "real(sp/dp) :: correlation" Correlation between \f$x\f$ and \f$y\f$.
227 : !> \param[in] "logical, optional :: mask(:)" 1D-array of logical values with `size(x)`.
228 : !! If present, only those locations in dat
229 : !! corresponding to the true values in mask are used.
230 :
231 : !> \note Input values must be floating points.
232 :
233 : !> \author Matthias Cuntz
234 : !> \date Nov 2011
235 : !> \date Dec 2011
236 : !! - covariance as <(x-<x>)(y-<y>)> instead of <xy>-<x><y>
237 :
238 : ! ------------------------------------------------------------------
239 :
240 : INTERFACE correlation
241 : MODULE PROCEDURE correlation_sp, correlation_dp
242 : END INTERFACE correlation
243 :
244 : ! ------------------------------------------------------------------
245 :
246 : !> \brief Covariance between vectors
247 :
248 : !> \details
249 : !! Calculates the covariance between two input vectors:\n
250 : !!
251 : !! \f[Cov(x,y) = \sum_i\frac{(x_i-\bar x)(y_i-\bar y)}{N}\f]
252 : !!
253 : !! If an optional mask is given, the average is only over those locations that correspond to true values in the mask.
254 : !! \f$x\f$ can be single or double precision. The result will have the same numerical precision.
255 : !!
256 : !! \b Example
257 : !!
258 : !! vec1 = (/ 1., 2, 3., -999., 5., 6. /)
259 : !! vec2 = (/ 1.3, 2.7, 3.9, 5.1, 6., 8. /)
260 : !! m = covariance(vec1, vec2, mask=((vec1 >= 0.) .and. (vec2 >= 0.)))
261 : !!
262 : !! See also example in test directory.
263 :
264 : !> \param[in] "real(sp/dp) :: x(:)" First 1D-array with input numbers.
265 : !> \param[in] "real(sp/dp) :: y(:)" Second 1D-array with input numbers.
266 : !> \param[in] "logical, optional :: mask(:)" 1D-array of logical values with size(x).
267 : !! If present, only those locations in dat
268 : !! corresponding to the true values in mask are used.
269 : !> \retval "real(sp/dp) :: covariance" Covariance between x and y.
270 :
271 : !> \note Input values must be floating points.
272 :
273 : !> \author Matthias Cuntz
274 : !> \date Nov 2011
275 : !> \date Dec 2011
276 : !! - covariance as <(x-<x>)(y-<y>)> instead of <xy>-<x><y>
277 :
278 : ! ------------------------------------------------------------------
279 :
280 : INTERFACE covariance
281 : MODULE PROCEDURE covariance_sp, covariance_dp
282 : END INTERFACE covariance
283 :
284 : ! ------------------------------------------------------------------
285 :
286 : !> \brief Kurtosis of a vector.
287 :
288 : !> \details
289 : !! Calculates the kurtosis of a vector, also called excess kurtosis:
290 : !!
291 : !! \f[ Kurt(x) = \verb|central_moment(x,4) / variance(x)**2 - 3|
292 : !! = \sum_i\frac{1}{N}\left(\frac{(x_i-\bar x)^2}{\mu_2(x)}\right)^2 - 3\f]
293 : !!
294 : !! If an optional mask is given, the average is only over those locations that correspond to true values in the mask.
295 : !! \f$x\f$ can be single or double precision. The result will have the same numerical precision.
296 : !!
297 : !! \b Example
298 : !!
299 : !! vec = (/ 1., 2, 3., -999., 5., 6. /)
300 : !! m = kurtosis(vec, mask=(vec >= 0.))
301 : !!
302 : !! See also example in test directory.
303 :
304 :
305 : !> \param[in] "real(sp/dp) :: dat(:)" 1D-array with input numbers.
306 : !> \param[in] "logical, optional :: mask(:)" 1D-array of logical values with `size(dat)`.
307 : !! If present, only those locations in dat
308 : !! corresponding to the true values in mask are used.
309 : !> \retval "real(sp/dp) :: kurtosis" Kurtosis of all elements in dat.
310 :
311 : !> \note Input values must be floating points.
312 :
313 : !> \authors Matthias Cuntz
314 : !> \date Nov 2011
315 :
316 : ! ------------------------------------------------------------------
317 :
318 : INTERFACE kurtosis
319 : MODULE PROCEDURE kurtosis_sp, kurtosis_dp
320 : END INTERFACE kurtosis
321 :
322 : ! ------------------------------------------------------------------
323 :
324 : !> \brief Mean of a vector.
325 :
326 : !> \details
327 : !! Calculates the average value of a vector, i.e. the first moment of a series of numbers:
328 : !!
329 : !! \f[\bar x = sum_i \frac{x_i}{N}
330 : !!
331 : !! If an optional mask is given, the mean is only over those locations that correspond to true values in the mask.
332 : !! \f$x\f$ can be single or double precision. The result will have the same numerical precision.
333 : !!
334 : !! \b Example
335 : !! \code{.f90}
336 : !! vec = (/ 1., 2, 3., -999., 5., 6. /)
337 : !! m = mean(vec, mask=(vec >= 0.))
338 : !! \endcode
339 : !! See also example in test directory.
340 : !!
341 : !> \param[in] "real(sp/dp) :: dat(:)" 1D-array with input numbers.
342 : !> \param[in] "logical, optional :: mask(:)" 1D-array of logical values with `size(dat)`.
343 : !! If present, only those locations in dat
344 : !! corresponding to the true values in mask are used.
345 : !> \retval "real(sp/dp) :: mean" Average of all elements in dat.
346 :
347 : !> \note Input values must be floating points.
348 :
349 : !> \author Matthias Cuntz
350 : !> \date Nov 2011
351 :
352 : ! ------------------------------------------------------------------
353 :
354 : INTERFACE mean
355 : MODULE PROCEDURE mean_sp, mean_dp
356 : END INTERFACE mean
357 :
358 : ! ------------------------------------------------------------------
359 :
360 : !> \brief R-s mixed central moment between vectors.
361 :
362 : !> \details
363 : !! Calculates the r,s-th mixed central moment between two vectors:
364 : !!
365 : !! \f[\mu_{r-s} = \sum_i\frac{(x_i-\bar{x})^r(y_i-\bar{y})^s}{N}\f]
366 : !!
367 : !! Note that covariace(x,y) = mixed_central_moment(x,y,1,1)
368 : !!
369 : !! If an optinal mask is given, the average is only over those locations that correspond to true values in the mask.
370 : !! \f$x\f$ can be single or double precision. The result will have the same numerical precision.
371 : !!
372 : !! \b Example
373 : !!
374 : !! vec1 = (/ 1., 2, 3., -999., 5., 6. /)
375 : !! vec2 = (/ 1.3, 2.7, 3.9, 5.1, 6., 8. /)
376 : !! m = mixed_central_moment(vec1, vec2, 2, 2, mask=((vec1 >= 0.) .and. (vec2 >= 0.)))
377 : !!
378 : !! See also example in test directory.
379 : !!
380 : !! \b Literature
381 : !!
382 : !! 1. LH Benedict & RD Gould, _Towards better uncertainty estimates for turbulence statistics_,
383 : !! Experiments in Fluids 22, 129-136, 1996
384 : !!
385 : !> \param[in] "real(sp/dp) :: x(:)" First 1D-array
386 : !> \param[in] "real(sp/dp) :: y(:)" Second 1D-array
387 : !> \param[in] "integer(i4) :: r" Order of x
388 : !> \param[in] "integer(i4) :: s" Order of y
389 : !> \param[in] "logical, optional :: mask(:)" 1D-array of logical values with size(x).
390 : !! If present, only those locations in dat
391 : !! corresponding to the true values in mask are used.
392 : !> \retval "real(sp/dp) :: mixed_central_moment" r,s-th mixed central moment between x and y
393 :
394 : !> \note Input values must be floating points.
395 :
396 : !> \author Matthias Cuntz
397 : !> \date Nov 2011
398 :
399 : ! ------------------------------------------------------------------
400 :
401 : INTERFACE mixed_central_moment
402 : MODULE PROCEDURE mixed_central_moment_sp, mixed_central_moment_dp
403 : END INTERFACE mixed_central_moment
404 :
405 : ! ------------------------------------------------------------------
406 :
407 : !> \brief Mixed central moment variance.
408 :
409 : !> \details
410 : !! Calculates the sample variance of r,s-th mixed central moment between two vectors.
411 : !! `mixed_central_moment_var` is something like the "error variance" of
412 : !! the r,s-th order mixed central moment sampling statistic.
413 : !!
414 : !! If an optinal mask is given, the average is only over those locations that correspond to true values in the mask.
415 : !! \f$x\f$ can be single or double precision. The result will have the same numerical precision.
416 : !!
417 : !! \b Example
418 : !!
419 : !! vec1 = (/ 1., 2, 3., -999., 5., 6. /)
420 : !! vec2 = (/ 1.3, 2.7, 3.9, 5.1, 6., 8. /)
421 : !! em = mixed_central_moment_var(vec1, vec2, 2, 2, mask=((vec1 >= 0.) .and. (vec2 >= 0.)))
422 : !!
423 : !! See also example in test directory.
424 : !!
425 : !! \b Literature
426 : !!
427 : !! 1. LH Benedict & RD Gould, _Towards better uncertainty estimates for turbulence statistics_,
428 : !! Experiments in Fluids 22, 129-136, 1996
429 : !!
430 : !> \param[in] "real(sp/dp) :: x(:)" First 1D-array
431 : !> \param[in] "real(sp/dp) :: y(:)" Second 1D-array
432 : !> \param[in] "integer(i4) :: r" Order of x
433 : !> \param[in] "integer(i4) :: s" Order of y
434 : !> \param[in] "logical, optional :: mask(:)" 1D-array of logical values with size(x).
435 : !! If present, only those locations in dat
436 : !! corresponding to the true values in mask are used.
437 : !> \retval "real(sp/dp) :: mixed_central_moment_var" Sampling variance of r,s-th mixed central moment between x and y
438 :
439 : !> \note Input values must be floating points.
440 :
441 : !> \author Matthias Cuntz
442 : !> \date Nov 2011
443 :
444 : INTERFACE mixed_central_moment_var
445 : MODULE PROCEDURE mixed_central_moment_var_sp, mixed_central_moment_var_dp
446 : END INTERFACE mixed_central_moment_var
447 :
448 : ! ------------------------------------------------------------------
449 :
450 : !> \brief First four moments, stddev and mean absolute devation.
451 :
452 : !> \details
453 : !! Calculates the first four sample/population moments of a vector, i.e. mean, variance, skewness, and kurtosis,
454 : !! as well as standard deviation and mean absolute devation.\n
455 : !!
456 : !! All output is optional.
457 : !!
458 : !! If an optinal mask is given, the calculations are over those locations that correspond to true values in the mask.
459 : !! \f$x\f$ can be single or double precision. The result will have the same numerical precision.
460 : !!
461 : !! \b Example
462 : !!
463 : !! vec = (/ 1., 2, 3., 4., 5., 6. /)
464 : !! call moment(vec, average=average, variance=variance, skewness=skewness, kurtosis=kurtosis, &
465 : !! mean=mean, stddev=stddev, absdev=absdev, mask=mask, sample=sample)
466 : !! --> average = 3.5
467 :
468 : !> \param[in] "real(sp/dp) :: dat(:)" 1D-array with input numbers.
469 : !> \param[in] "logical, optional :: mask(:)" 1D-array of logical values with `size(dat)`.
470 : !! If present, only those locations in vec
471 : !! corresponding to the true values in mask are used.
472 : !> \param[in] "logical, optional :: sample" Logical value.
473 : !! If present and False, the population variance and
474 : !! std-dev will be calculated (divide by n).
475 :
476 : !> \param[out] "real(sp/dp), optional :: average" Average of input vector.
477 : !> \param[out] "real(sp/dp), optional :: variance" Sample variance of input vector (either a sample or pupulation moment).
478 : !> \param[out] "real(sp/dp), optional :: skewness" Skewness of input vector.
479 : !> \param[out] "real(sp/dp), optional :: kurtosis" Excess kurtosis of input vector.
480 : !> \param[out] "real(sp/dp), optional :: mean" Same as average.
481 : !> \param[out] "real(sp/dp), optional :: stddev" Sqaure root of variance (either a sample or pupulation moment).
482 : !> \param[out] "real(sp/dp), optional :: absdev" Mean absolute deviations from average.
483 :
484 : !> \note Input values must be floating points. Inpt and all optional outputs must have same kind.
485 :
486 : !> \author Matthias Cuntz
487 : !> \date Nov 2011
488 : !> \author Sebastian Mueller
489 : !> \date Dec 2019
490 : !! - added optional sample input-para to switch sample to population variance and std-dev.
491 :
492 : ! ------------------------------------------------------------------
493 :
494 : INTERFACE moment
495 : MODULE PROCEDURE moment_sp, moment_dp
496 : END INTERFACE moment
497 :
498 : ! ------------------------------------------------------------------
499 :
500 : !> \brief Skewness of a vector
501 :
502 : !> \details
503 : !! Calculates the skewness of a vector, i.e. the third standardised moment:
504 : !!
505 : !! \f[\tilde \mu_3 = \sum_i\left(\frac{(x-\bar x)}{\sigma_x}\right)^3\f]
506 : !!
507 : !! If an optinal mask is given, the average is only over those locations that correspond to true values in the mask.
508 : !! \f$x\f$ can be single or double precision. The result will have the same numerical precision.
509 : !!
510 : !! \b Example
511 : !!
512 : !! vec = (/ 1., 2, 3., -999., 5., 6. /)
513 : !! m = skewness(vec, mask=(vec >= 0.))
514 : !!
515 : !! See also example in test directory.
516 :
517 :
518 : !> \param[in] "real(sp/dp) :: dat(:)" 1D-array with input numbers.
519 : !> \param[in] "logical, optional :: mask(:)" 1D-array of logical values with `size(dat)`.
520 : !! If present, only those locations in vec
521 : !! corresponding to the true values in mask are used.
522 : !> \retval "real(sp/dp) :: skewness" Skewness of all elements in dat.
523 :
524 : !> \note Input values must be floating points.
525 :
526 : !> \author Matthias Cuntz
527 : !> \date Nov 2011
528 : INTERFACE skewness
529 : MODULE PROCEDURE skewness_sp, skewness_dp
530 : END INTERFACE skewness
531 :
532 : ! ------------------------------------------------------------------
533 :
534 : !> \brief Standard deviation of a vector.
535 :
536 : !> \details
537 : !! Calculates the sample standard deviation of a vector, i.e. the square root of the second moment
538 : !!
539 : !! \f[\sigma_x = \sqrt{\sum_i\frac{(x_i-\bar x)^2}{N-1}}\f]
540 : !!
541 : !! If an optinal mask is given, the average is only over those locations that correspond to true values in the mask.
542 : !! \f$x\f$ can be single or double precision. The result will have the same numerical precision.
543 : !!
544 : !! \b Example
545 : !!
546 : !! vec = (/ 1., 2, 3., -999., 5., 6. /)
547 : !! m = stddev(vec, mask=(vec >= 0.))
548 : !!
549 : !! See also example in test directory
550 :
551 : !> \param[in] "real(sp/dp) :: dat(:)" 1D-array with input numbers.
552 : !> \param[in] "logical, optional :: mask(:)" 1D-array of logical values with `size(dat)`.
553 : !! If present, only those locations in vec
554 : !! corresponding to the true values in mask are used.
555 : !> \retval "real(sp/dp) :: stddev" Sample standard deviation of all elements in dat.
556 :
557 : !> \note
558 : !! Input values must be floating points.\n
559 : !! This is the sample standard deviation. The population standard deviation is:
560 : !! `popstddev = stddev * sqrt((n-1)/n)`
561 :
562 : !> \author Matthias Cuntz
563 : !> \date Nov 2011
564 :
565 : ! ------------------------------------------------------------------
566 :
567 : INTERFACE stddev
568 : MODULE PROCEDURE stddev_sp, stddev_dp
569 : END INTERFACE stddev
570 :
571 : ! ------------------------------------------------------------------
572 :
573 : !> \brief Standard deviation of a vector.
574 :
575 : !> \details
576 : !! Calculates the sample variance of a vector, i.e. the second moment
577 : !!
578 : !! \f[\sigma_x^2 = \sum_i\frac{(x_i-\bar x)^2}{N-1}\f]
579 : !!
580 : !! If an optional mask is given, the average is only over those locations that correspond to true values in the mask.
581 : !! \f$x\f$ can be single or double precision. The result will have the same numerical precision.
582 :
583 : !!
584 : !! \b Example
585 : !!
586 : !! vec = (/ 1., 2, 3., -999., 5., 6. /)
587 : !! m = variance(vec, mask=(vec >= 0.))
588 : !!
589 : !! See also example in test directory
590 :
591 : !> \param[in] "real(sp/dp) :: dat(:)" 1D-array with input numbers.
592 : !> \param[in] "logical, optional :: mask(:)" 1D-array of logical values with `size(dat)`.
593 : !! If present, only those locations in vec
594 : !! corresponding to the true values in mask are used.
595 : !> \retval "real(sp/dp) :: variance" Sample variance of all elements in dat.
596 :
597 : !> \note
598 : !! Input values must be floating points.\n
599 : !! This is the sample variance. The population variance is:
600 : !! `var = variance * (n-1)/n`
601 :
602 : !> \author Matthias Cuntz
603 : !> \date Nov 2011
604 : INTERFACE variance
605 : MODULE PROCEDURE variance_sp, variance_dp
606 : END INTERFACE variance
607 :
608 : ! ------------------------------------------------------------------
609 :
610 : PRIVATE
611 :
612 : ! ------------------------------------------------------------------
613 :
614 : CONTAINS
615 :
616 : ! ------------------------------------------------------------------
617 :
618 6 : FUNCTION absdev_dp(dat, mask)
619 :
620 : IMPLICIT NONE
621 :
622 : REAL(dp), DIMENSION(:), INTENT(IN) :: dat
623 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
624 : REAL(dp) :: absdev_dp
625 :
626 3 : REAL(dp) :: n
627 :
628 3 : REAL(dp) :: ave
629 3 : LOGICAL, DIMENSION(size(dat)) :: maske
630 :
631 3 : if (present(mask)) then
632 1 : if (size(mask) .ne. size(dat)) stop 'Error absdev_dp: size(mask) .ne. size(dat)'
633 12 : maske = mask
634 11 : n = real(count(maske), dp)
635 : else
636 21 : maske(:) = .true.
637 2 : n = real(size(dat), dp)
638 : end if
639 3 : if (n .le. (1.0_dp + tiny(1.0_dp))) stop 'absdev_dp: n must be at least 2'
640 :
641 : ! Average
642 35 : ave = sum(dat(:), mask = maske) / n
643 : ! Sum of absolute deviation
644 32 : absdev_dp = sum(abs(dat(:) - ave), mask = maske) / n
645 :
646 3 : END FUNCTION absdev_dp
647 :
648 :
649 6 : FUNCTION absdev_sp(dat, mask)
650 :
651 : IMPLICIT NONE
652 :
653 : REAL(sp), DIMENSION(:), INTENT(IN) :: dat
654 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
655 : REAL(sp) :: absdev_sp
656 :
657 3 : REAL(sp) :: n
658 :
659 3 : REAL(sp) :: ave
660 3 : LOGICAL, DIMENSION(size(dat)) :: maske
661 :
662 3 : if (present(mask)) then
663 1 : if (size(mask) .ne. size(dat)) stop 'Error absdev_sp: size(mask) .ne. size(dat)'
664 12 : maske = mask
665 11 : n = real(count(maske), sp)
666 : else
667 21 : maske(:) = .true.
668 2 : n = real(size(dat), sp)
669 : end if
670 3 : if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'absdev_sp: n must be at least 2'
671 :
672 : ! Average
673 35 : ave = sum(dat(:), mask = maske) / n
674 : ! Sum of absolute deviation
675 32 : absdev_sp = sum(abs(dat(:) - ave), mask = maske) / n
676 :
677 3 : END FUNCTION absdev_sp
678 :
679 : ! ------------------------------------------------------------------
680 :
681 130 : FUNCTION average_dp(dat, mask)
682 :
683 : IMPLICIT NONE
684 :
685 : REAL(dp), DIMENSION(:), INTENT(IN) :: dat
686 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
687 : REAL(dp) :: average_dp
688 :
689 65 : REAL(dp) :: n
690 65 : LOGICAL, DIMENSION(size(dat)) :: maske
691 :
692 65 : if (present(mask)) then
693 59 : if (size(mask) .ne. size(dat)) stop 'Error average_dp: size(mask) .ne. size(dat)'
694 8948 : maske = mask
695 8889 : n = real(count(maske), dp)
696 : else
697 762 : maske(:) = .true.
698 6 : n = real(size(dat), dp)
699 : end if
700 :
701 : ! Average
702 9716 : average_dp = sum(dat(:), mask = maske) / n
703 :
704 3 : END FUNCTION average_dp
705 :
706 :
707 130 : FUNCTION average_sp(dat, mask)
708 :
709 : IMPLICIT NONE
710 :
711 : REAL(sp), DIMENSION(:), INTENT(IN) :: dat
712 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
713 : REAL(sp) :: average_sp
714 :
715 65 : REAL(sp) :: n
716 65 : LOGICAL, DIMENSION(size(dat)) :: maske
717 :
718 65 : if (present(mask)) then
719 59 : if (size(mask) .ne. size(dat)) stop 'Error average_sp: size(mask) .ne. size(dat)'
720 8948 : maske = mask
721 8889 : n = real(count(maske), sp)
722 : else
723 762 : maske(:) = .true.
724 6 : n = real(size(dat), sp)
725 : end if
726 :
727 : ! Average
728 9716 : average_sp = sum(dat(:), mask = maske) / n
729 :
730 65 : END FUNCTION average_sp
731 :
732 : ! ------------------------------------------------------------------
733 :
734 36 : FUNCTION central_moment_dp(x, r, mask)
735 :
736 : IMPLICIT NONE
737 :
738 : REAL(dp), DIMENSION(:), INTENT(IN) :: x
739 : INTEGER(i4), INTENT(IN) :: r
740 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
741 : REAL(dp) :: central_moment_dp
742 :
743 18 : REAL(dp) :: n, mx
744 18 : LOGICAL, DIMENSION(size(x)) :: maske
745 :
746 18 : if (r .lt. 0) then
747 18 : central_moment_dp = 0.0_dp
748 : return
749 : end if
750 18 : if (r .eq. 0) then
751 18 : central_moment_dp = 1.0_dp
752 : return
753 : end if
754 :
755 18 : if (present(mask)) then
756 16 : if (size(mask) .ne. size(x)) stop 'Error central_moment_dp: size(mask) .ne. size(x)'
757 171 : maske = mask
758 171 : n = real(count(maske), dp)
759 : else
760 21 : maske(:) = .true.
761 2 : n = real(size(x), dp)
762 : end if
763 18 : if (n .le. (2.0_dp + tiny(2.0_dp))) stop 'central_moment_dp: n must be at least 3'
764 :
765 : ! average
766 192 : mx = sum(x, mask = maske) / n
767 : ! central moment
768 192 : central_moment_dp = sum((x - mx)**r, mask = maske) / n
769 :
770 83 : END FUNCTION central_moment_dp
771 :
772 :
773 36 : FUNCTION central_moment_sp(x, r, mask)
774 :
775 : IMPLICIT NONE
776 :
777 : REAL(sp), DIMENSION(:), INTENT(IN) :: x
778 : INTEGER(i4), INTENT(IN) :: r
779 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
780 : REAL(sp) :: central_moment_sp
781 :
782 18 : REAL(sp) :: n, mx
783 18 : LOGICAL, DIMENSION(size(x)) :: maske
784 :
785 18 : if (r .lt. 0) then
786 18 : central_moment_sp = 0.0_sp
787 : return
788 : end if
789 18 : if (r .eq. 0) then
790 18 : central_moment_sp = 1.0_sp
791 : return
792 : end if
793 :
794 18 : if (present(mask)) then
795 16 : if (size(mask) .ne. size(x)) stop 'Error central_moment_sp: size(mask) .ne. size(x)'
796 171 : maske = mask
797 171 : n = real(count(maske), sp)
798 : else
799 21 : maske(:) = .true.
800 2 : n = real(size(x), sp)
801 : end if
802 18 : if (n .le. (2.0_sp + tiny(2.0_sp))) stop 'central_moment_sp: n must be at least 3'
803 :
804 : ! average
805 192 : mx = sum(x, mask = maske) / n
806 : ! central moment
807 192 : central_moment_sp = sum((x - mx)**r, mask = maske) / n
808 :
809 36 : END FUNCTION central_moment_sp
810 :
811 : ! ------------------------------------------------------------------
812 :
813 6 : FUNCTION central_moment_var_dp(x, r, mask)
814 :
815 : IMPLICIT NONE
816 :
817 : REAL(dp), DIMENSION(:), INTENT(IN) :: x
818 : INTEGER(i4), INTENT(IN) :: r
819 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
820 : REAL(dp) :: central_moment_var_dp
821 :
822 3 : REAL(dp) :: n, rr, u2r, ur, urm1, urp1, u2
823 3 : LOGICAL, DIMENSION(size(x)) :: maske
824 :
825 3 : if (r.le.1) then
826 3 : central_moment_var_dp = 0.0_dp
827 : return
828 : end if
829 :
830 3 : if (present(mask)) then
831 1 : if (size(mask) .ne. size(x)) stop 'Error central_moment_var_dp: size(mask) .ne. size(x)'
832 11 : maske = mask
833 11 : n = real(count(maske), dp)
834 : else
835 21 : maske(:) = .true.
836 2 : n = real(size(x), dp)
837 : end if
838 3 : if (n .le. (2.0_dp + tiny(2.0_dp))) stop 'central_moment_var_dp: n must be at least 3'
839 :
840 3 : u2r = central_moment(x, 2 * r, mask = maske)
841 3 : ur = central_moment(x, r, mask = maske)
842 3 : urm1 = central_moment(x, r - 1, mask = maske)
843 3 : u2 = central_moment(x, 2, mask = maske)
844 3 : urp1 = central_moment(x, r + 1, mask = maske)
845 3 : rr = real(r, dp)
846 3 : central_moment_var_dp = (u2r - ur * ur + rr * rr * urm1 * urm1 * u2 - 2.0_dp * rr * urp1 * urm1) / n
847 :
848 21 : END FUNCTION central_moment_var_dp
849 :
850 :
851 6 : FUNCTION central_moment_var_sp(x, r, mask)
852 :
853 : IMPLICIT NONE
854 :
855 : REAL(sp), DIMENSION(:), INTENT(IN) :: x
856 : INTEGER(i4), INTENT(IN) :: r
857 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
858 : REAL(sp) :: central_moment_var_sp
859 :
860 3 : REAL(sp) :: n, rr, u2r, ur, urm1, urp1, u2
861 3 : LOGICAL, DIMENSION(size(x)) :: maske
862 :
863 3 : if (r.le.1) then
864 3 : central_moment_var_sp = 0.0_sp
865 : return
866 : end if
867 :
868 3 : if (present(mask)) then
869 1 : if (size(mask) .ne. size(x)) stop 'Error central_moment_var_sp: size(mask) .ne. size(x)'
870 11 : maske = mask
871 11 : n = real(count(maske), sp)
872 : else
873 21 : maske(:) = .true.
874 2 : n = real(size(x), sp)
875 : end if
876 3 : if (n .le. (2.0_sp + tiny(2.0_sp))) stop 'central_moment_var_sp: n must be at least 3'
877 :
878 3 : u2r = central_moment(x, 2 * r, mask = maske)
879 3 : ur = central_moment(x, r, mask = maske)
880 3 : urm1 = central_moment(x, r - 1, mask = maske)
881 3 : u2 = central_moment(x, 2, mask = maske)
882 3 : urp1 = central_moment(x, r + 1, mask = maske)
883 3 : rr = real(r, sp)
884 3 : central_moment_var_sp = (u2r - ur * ur + rr * rr * urm1 * urm1 * u2 - 2.0_sp * rr * urp1 * urm1) / n
885 :
886 6 : END FUNCTION central_moment_var_sp
887 :
888 : ! ------------------------------------------------------------------
889 :
890 18 : FUNCTION correlation_dp(x, y, mask)
891 :
892 : IMPLICIT NONE
893 :
894 : REAL(dp), DIMENSION(:), INTENT(IN) :: x
895 : REAL(dp), DIMENSION(:), INTENT(IN) :: y
896 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
897 : REAL(dp) :: correlation_dp
898 :
899 9 : REAL(dp) :: n
900 9 : REAL(dp) :: mx, my
901 9 : REAL(dp) :: sx, sy, covar
902 9 : LOGICAL, DIMENSION(size(x)) :: maske
903 :
904 9 : if (size(x) .ne. size(y)) stop 'Error correlation_dp: size(x) .ne. size(y)'
905 9 : if (present(mask)) then
906 7 : if (size(mask) .ne. size(x)) stop 'Error correlation_dp: size(mask) .ne. size(x)'
907 1104 : maske = mask
908 1097 : n = real(count(maske), dp)
909 : else
910 21 : maske(:) = .true.
911 2 : n = real(size(x), dp)
912 : end if
913 9 : if (n .le. (1.0_dp + tiny(1.0_dp))) stop 'correlation_dp: n must be at least 2'
914 :
915 : ! Mean and Stddev of x and y
916 9 : call moment(x, mx, stddev = sx, mask = maske)
917 9 : call moment(y, my, stddev = sy, mask = maske)
918 1118 : covar = sum((x - mx) * (y - my), mask = maske) / n
919 : ! correlation
920 9 : correlation_dp = covar / (sx * sy)
921 :
922 3 : END FUNCTION correlation_dp
923 :
924 :
925 18 : FUNCTION correlation_sp(x, y, mask)
926 :
927 : IMPLICIT NONE
928 :
929 : REAL(sp), DIMENSION(:), INTENT(IN) :: x
930 : REAL(sp), DIMENSION(:), INTENT(IN) :: y
931 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
932 : REAL(sp) :: correlation_sp
933 :
934 9 : REAL(sp) :: n
935 9 : REAL(sp) :: mx, my
936 9 : REAL(sp) :: sx, sy, covar
937 9 : LOGICAL, DIMENSION(size(x)) :: maske
938 :
939 9 : if (size(x) .ne. size(y)) stop 'Error correlation_sp: size(x) .ne. size(y)'
940 9 : if (present(mask)) then
941 7 : if (size(mask) .ne. size(x)) stop 'Error correlation_sp: size(mask) .ne. size(x)'
942 1104 : maske = mask
943 1097 : n = real(count(maske), sp)
944 : else
945 21 : maske(:) = .true.
946 2 : n = real(size(x), sp)
947 : end if
948 9 : if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'correlation_sp: n must be at least 2'
949 :
950 : ! Mean and Stddev of x and y
951 9 : call moment(x, mx, stddev = sx, mask = maske)
952 9 : call moment(y, my, stddev = sy, mask = maske)
953 1118 : covar = sum((x - mx) * (y - my), mask = maske) / n
954 : ! correlation
955 9 : correlation_sp = covar / (sx * sy)
956 :
957 9 : END FUNCTION correlation_sp
958 :
959 : ! ------------------------------------------------------------------
960 :
961 6 : FUNCTION covariance_dp(x, y, mask)
962 :
963 : IMPLICIT NONE
964 :
965 : REAL(dp), DIMENSION(:), INTENT(IN) :: x
966 : REAL(dp), DIMENSION(:), INTENT(IN) :: y
967 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
968 : REAL(dp) :: covariance_dp
969 :
970 3 : REAL(dp) :: n
971 3 : REAL(dp) :: mx, my
972 3 : LOGICAL, DIMENSION(size(x)) :: maske
973 :
974 3 : if (size(x) .ne. size(y)) stop 'Error covariance_dp: size(x) .ne. size(y)'
975 3 : if (present(mask)) then
976 1 : if (size(mask) .ne. size(x)) stop 'Error covariance_dp: size(mask) .ne. size(x)'
977 12 : maske = mask
978 11 : n = real(count(maske), dp)
979 : else
980 21 : maske(:) = .true.
981 2 : n = real(size(x), dp)
982 : end if
983 3 : if (n .le. (1.0_dp + tiny(1.0_dp))) stop 'covariance_dp: n must be at least 2'
984 :
985 : ! Mean of x and y
986 3 : mx = mean(x, mask = maske)
987 3 : my = mean(y, mask = maske)
988 35 : covariance_dp = sum((x - mx) * (y - my), mask = maske) / n
989 :
990 9 : END FUNCTION covariance_dp
991 :
992 :
993 6 : FUNCTION covariance_sp(x, y, mask)
994 :
995 : IMPLICIT NONE
996 :
997 : REAL(sp), DIMENSION(:), INTENT(IN) :: x
998 : REAL(sp), DIMENSION(:), INTENT(IN) :: y
999 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
1000 : REAL(sp) :: covariance_sp
1001 :
1002 3 : REAL(sp) :: n
1003 3 : REAL(sp) :: mx, my
1004 3 : LOGICAL, DIMENSION(size(x)) :: maske
1005 :
1006 3 : if (size(x) .ne. size(y)) stop 'Error covariance_sp: size(x) .ne. size(y)'
1007 3 : if (present(mask)) then
1008 1 : if (size(mask) .ne. size(x)) stop 'Error covariance_sp: size(mask) .ne. size(x)'
1009 12 : maske = mask
1010 11 : n = real(count(maske), sp)
1011 : else
1012 21 : maske(:) = .true.
1013 2 : n = real(size(x), sp)
1014 : end if
1015 3 : if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'covariance_sp: n must be at least 2'
1016 :
1017 : ! Mean of x and y
1018 3 : mx = mean(x, mask = maske)
1019 3 : my = mean(y, mask = maske)
1020 35 : covariance_sp = sum((x - mx) * (y - my), mask = maske) / n
1021 :
1022 3 : END FUNCTION covariance_sp
1023 :
1024 : ! ------------------------------------------------------------------
1025 :
1026 6 : FUNCTION kurtosis_dp(dat, mask)
1027 :
1028 : IMPLICIT NONE
1029 :
1030 : REAL(dp), DIMENSION(:), INTENT(IN) :: dat
1031 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
1032 : REAL(dp) :: kurtosis_dp
1033 :
1034 3 : REAL(dp) :: n
1035 :
1036 35 : REAL(dp) :: ep, ave, var
1037 67 : REAL(dp), DIMENSION(size(dat)) :: p, s
1038 3 : LOGICAL, DIMENSION(size(dat)) :: maske
1039 :
1040 3 : if (present(mask)) then
1041 1 : if (size(mask) .ne. size(dat)) stop 'Error kurtosis_dp: size(mask) .ne. size(dat)'
1042 12 : maske = mask
1043 11 : n = real(count(maske), dp)
1044 : else
1045 21 : maske(:) = .true.
1046 2 : n = real(size(dat), dp)
1047 : end if
1048 3 : if (n .le. (1.0_dp + tiny(1.0_dp))) stop 'kurtosis_dp: n must be at least 2'
1049 :
1050 : ! Average
1051 35 : ave = sum(dat(:), mask = maske) / n
1052 32 : s(:) = dat(:) - ave
1053 : ! Variance / Standard deviation
1054 32 : ep = sum(s(:), mask = maske)
1055 32 : p(:) = s(:) * s(:)
1056 32 : var = sum(p(:), mask = maske)
1057 3 : var = (var - ep * ep / n) / (n - 1.0_dp)
1058 3 : if (abs(var) .lt. tiny(0.0_dp)) stop 'kurtosis_dp: no kurtosis when zero variance'
1059 : ! Kurtosis
1060 32 : p(:) = p(:) * s(:) * s(:)
1061 32 : kurtosis_dp = sum(p(:), mask = maske)
1062 3 : kurtosis_dp = kurtosis_dp / (n * var * var) - 3.0_dp
1063 :
1064 3 : END FUNCTION kurtosis_dp
1065 :
1066 :
1067 6 : FUNCTION kurtosis_sp(dat, mask)
1068 :
1069 : IMPLICIT NONE
1070 :
1071 : REAL(sp), DIMENSION(:), INTENT(IN) :: dat
1072 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
1073 : REAL(sp) :: kurtosis_sp
1074 :
1075 3 : REAL(sp) :: n
1076 :
1077 35 : REAL(sp) :: ep, ave, var
1078 67 : REAL(sp), DIMENSION(size(dat)) :: p, s
1079 3 : LOGICAL, DIMENSION(size(dat)) :: maske
1080 :
1081 3 : if (present(mask)) then
1082 1 : if (size(mask) .ne. size(dat)) stop 'Error kurtosis_sp: size(mask) .ne. size(dat)'
1083 12 : maske = mask
1084 11 : n = real(count(maske), sp)
1085 : else
1086 21 : maske(:) = .true.
1087 2 : n = real(size(dat), sp)
1088 : end if
1089 3 : if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'kurtosis_sp: n must be at least 2'
1090 :
1091 : ! Average
1092 35 : ave = sum(dat(:), mask = maske) / n
1093 32 : s(:) = dat(:) - ave
1094 : ! Variance / Standard deviation
1095 32 : ep = sum(s(:), mask = maske)
1096 32 : p(:) = s(:) * s(:)
1097 32 : var = sum(p(:), mask = maske)
1098 3 : var = (var - ep * ep / n) / (n - 1.0_sp)
1099 3 : if (abs(var) .lt. tiny(0.0_sp)) stop 'kurtosis_sp: no kurtosis when zero variance'
1100 : ! Kurtosis
1101 32 : p(:) = p(:) * s(:) * s(:)
1102 32 : kurtosis_sp = sum(p(:), mask = maske)
1103 3 : kurtosis_sp = kurtosis_sp / (n * var * var) - 3.0_sp
1104 :
1105 3 : END FUNCTION kurtosis_sp
1106 :
1107 : ! ------------------------------------------------------------------
1108 :
1109 9084348 : FUNCTION mean_dp(dat, mask)
1110 :
1111 : IMPLICIT NONE
1112 :
1113 : REAL(dp), DIMENSION(:), INTENT(IN) :: dat
1114 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
1115 : REAL(dp) :: mean_dp
1116 :
1117 4542174 : REAL(dp) :: n
1118 :
1119 4542174 : LOGICAL, DIMENSION(size(dat)) :: maske
1120 :
1121 4542174 : if (present(mask)) then
1122 4542163 : if (size(mask) .ne. size(dat)) stop 'Error mean_dp: size(mask) .ne. size(dat)'
1123 46356358 : maske = mask
1124 41814195 : n = real(count(maske), dp)
1125 : else
1126 225030 : maske(:) = .true.
1127 11 : n = real(size(dat), dp)
1128 : end if
1129 :
1130 : ! Mean
1131 46581399 : mean_dp = sum(dat(:), mask = maske) / n
1132 :
1133 3 : END FUNCTION mean_dp
1134 :
1135 : !> \copydoc mean
1136 18 : FUNCTION mean_sp(dat, mask)
1137 :
1138 : IMPLICIT NONE
1139 :
1140 : REAL(sp), DIMENSION(:), INTENT(IN) :: dat
1141 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
1142 : REAL(sp) :: mean_sp
1143 :
1144 9 : REAL(sp) :: n
1145 :
1146 9 : LOGICAL, DIMENSION(size(dat)) :: maske
1147 :
1148 9 : if (present(mask)) then
1149 7 : if (size(mask) .ne. size(dat)) stop 'Error mean_sp: size(mask) .ne. size(dat)'
1150 82 : maske = mask
1151 75 : n = real(count(maske), sp)
1152 : else
1153 21 : maske(:) = .true.
1154 2 : n = real(size(dat), sp)
1155 : end if
1156 :
1157 : ! Mean
1158 105 : mean_sp = sum(dat(:), mask = maske) / n
1159 :
1160 4542174 : END FUNCTION mean_sp
1161 :
1162 : ! ------------------------------------------------------------------
1163 :
1164 60 : FUNCTION mixed_central_moment_dp(x, y, r, s, mask)
1165 :
1166 : IMPLICIT NONE
1167 :
1168 : REAL(dp), DIMENSION(:), INTENT(IN) :: x
1169 : REAL(dp), DIMENSION(:), INTENT(IN) :: y
1170 : INTEGER(i4), INTENT(IN) :: r
1171 : INTEGER(i4), INTENT(IN) :: s
1172 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
1173 : REAL(dp) :: mixed_central_moment_dp
1174 :
1175 30 : REAL(dp) :: n, mx, my
1176 610 : REAL(dp), DIMENSION(size(x)) :: xx, yy
1177 30 : LOGICAL, DIMENSION(size(x)) :: maske
1178 :
1179 30 : if (r.lt.0 .or. s.lt.0) then
1180 30 : mixed_central_moment_dp = 0.0_dp
1181 : return
1182 : end if
1183 :
1184 30 : if (size(x) .ne. size(y)) stop 'Error mixed_central_moment_dp: size(x) .ne. size(y)'
1185 30 : if (present(mask)) then
1186 28 : if (size(mask) .ne. size(x)) stop 'Error mixed_central_moment_dp: size(mask) .ne. size(x)'
1187 299 : maske = mask
1188 299 : n = real(count(maske), dp)
1189 : else
1190 21 : maske(:) = .true.
1191 2 : n = real(size(x), dp)
1192 : end if
1193 30 : if (n .le. (2.0_dp + tiny(2.0_dp))) stop 'mixed_central_moment_dp: n must be at least 3'
1194 :
1195 : ! Averages of x and y
1196 320 : mx = sum(x, mask = maske) / n
1197 320 : my = sum(y, mask = maske) / n
1198 : ! Mixed central moment
1199 30 : if (r>0) then
1200 256 : xx = (x - mx)**r
1201 : else
1202 64 : xx = 1._dp
1203 : end if
1204 30 : if (s>0) then
1205 256 : yy = (y - my)**s
1206 : else
1207 64 : yy = 1._dp
1208 : end if
1209 320 : mixed_central_moment_dp = sum(xx * yy, mask = maske) / n
1210 :
1211 39 : END FUNCTION mixed_central_moment_dp
1212 :
1213 :
1214 60 : FUNCTION mixed_central_moment_sp(x, y, r, s, mask)
1215 :
1216 : IMPLICIT NONE
1217 :
1218 : REAL(sp), DIMENSION(:), INTENT(IN) :: x
1219 : REAL(sp), DIMENSION(:), INTENT(IN) :: y
1220 : INTEGER(i4), INTENT(IN) :: r
1221 : INTEGER(i4), INTENT(IN) :: s
1222 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
1223 : REAL(sp) :: mixed_central_moment_sp
1224 :
1225 30 : REAL(sp) :: n, mx, my
1226 610 : REAL(sp), DIMENSION(size(x)) :: xx, yy
1227 30 : LOGICAL, DIMENSION(size(x)) :: maske
1228 :
1229 30 : if (r.lt.0 .or. s.lt.0) then
1230 30 : mixed_central_moment_sp = 0.0_sp
1231 : return
1232 : end if
1233 :
1234 30 : if (size(x) .ne. size(y)) stop 'Error mixed_central_moment_sp: size(x) .ne. size(y)'
1235 30 : if (present(mask)) then
1236 28 : if (size(mask) .ne. size(x)) stop 'Error mixed_central_moment_sp: size(mask) .ne. size(x)'
1237 299 : maske = mask
1238 299 : n = real(count(maske), sp)
1239 : else
1240 21 : maske(:) = .true.
1241 2 : n = real(size(x), sp)
1242 : end if
1243 30 : if (n .le. (2.0_sp + tiny(2.0_sp))) stop 'mixed_central_moment_sp: n must be at least 3'
1244 :
1245 : ! Averages of x and y
1246 320 : mx = sum(x, mask = maske) / n
1247 320 : my = sum(y, mask = maske) / n
1248 : ! Mixed central moment
1249 30 : if (r>0) then
1250 256 : xx = (x - mx)**r
1251 : else
1252 64 : xx = 1._sp
1253 : end if
1254 30 : if (s>0) then
1255 256 : yy = (y - my)**s
1256 : else
1257 64 : yy = 1._sp
1258 : end if
1259 320 : mixed_central_moment_sp = sum(xx * yy, mask = maske) / n
1260 :
1261 60 : END FUNCTION mixed_central_moment_sp
1262 :
1263 : ! ------------------------------------------------------------------
1264 :
1265 6 : FUNCTION mixed_central_moment_var_dp(x, y, r, s, mask)
1266 : ! Error variance of mixed central moment (Benedict & Gould 1996)
1267 : IMPLICIT NONE
1268 :
1269 : REAL(dp), DIMENSION(:), INTENT(IN) :: x
1270 : REAL(dp), DIMENSION(:), INTENT(IN) :: y
1271 : INTEGER(i4), INTENT(IN) :: r
1272 : INTEGER(i4), INTENT(IN) :: s
1273 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
1274 : REAL(dp) :: mixed_central_moment_var_dp
1275 :
1276 3 : REAL(dp) :: u2r2s, urs, urm1s, u20, urp1s, ursm1, u02, ursp1, u11
1277 3 : REAL(dp) :: n, rr, ss
1278 3 : LOGICAL, DIMENSION(size(x)) :: maske
1279 :
1280 3 : if (size(x) .ne. size(y)) stop 'Error mixed_central_moment_var_dp: size(x) .ne. size(y)'
1281 3 : if (present(mask)) then
1282 1 : if (size(mask) .ne. size(x)) stop 'Error mixed_central_moment_var_dp: size(mask) .ne. size(x)'
1283 11 : maske = mask
1284 11 : n = real(count(maske), dp)
1285 : else
1286 21 : maske(:) = .true.
1287 2 : n = real(size(x), dp)
1288 : end if
1289 3 : if (n .le. (2.0_dp + tiny(2.0_dp))) stop 'mixed_central_moment_var_dp: n must be at least 3'
1290 :
1291 3 : u2r2s = mixed_central_moment(x, y, 2 * r, 2 * s, mask = maske)
1292 3 : urs = mixed_central_moment(x, y, r, s, mask = maske)
1293 3 : urm1s = mixed_central_moment(x, y, r - 1, s, mask = maske)
1294 3 : u20 = mixed_central_moment(x, y, 2, 0, mask = maske)
1295 3 : urp1s = mixed_central_moment(x, y, r + 1, s, mask = maske)
1296 3 : ursm1 = mixed_central_moment(x, y, r, s - 1, mask = maske)
1297 3 : u02 = mixed_central_moment(x, y, 0, 2, mask = maske)
1298 3 : ursp1 = mixed_central_moment(x, y, r, s + 1, mask = maske)
1299 3 : u11 = mixed_central_moment(x, y, 1, 1, mask = maske)
1300 3 : rr = real(r, dp)
1301 3 : ss = real(s, dp)
1302 :
1303 : mixed_central_moment_var_dp = (u2r2s - urs * urs &
1304 : + rr * rr**u20 * urm1s * urm1s + ss * ss * u02 * ursm1 * ursm1 &
1305 : + 2.0_dp * rr * ss * u11 * urm1s * ursm1 &
1306 3 : - 2.0_dp * rr * urp1s * urm1s - 2.0_dp * ss * ursp1 * ursm1) / n
1307 :
1308 30 : END FUNCTION mixed_central_moment_var_dp
1309 :
1310 :
1311 6 : FUNCTION mixed_central_moment_var_sp(x, y, r, s, mask)
1312 : ! Error variance of mixed central moment (Benedict & Gould 1996)
1313 : IMPLICIT NONE
1314 :
1315 : REAL(sp), DIMENSION(:), INTENT(IN) :: x
1316 : REAL(sp), DIMENSION(:), INTENT(IN) :: y
1317 : INTEGER(i4), INTENT(IN) :: r
1318 : INTEGER(i4), INTENT(IN) :: s
1319 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
1320 : REAL(sp) :: mixed_central_moment_var_sp
1321 :
1322 3 : REAL(sp) :: u2r2s, urs, urm1s, u20, urp1s, ursm1, u02, ursp1, u11
1323 3 : REAL(sp) :: n, rr, ss
1324 3 : LOGICAL, DIMENSION(size(x)) :: maske
1325 :
1326 3 : if (size(x) .ne. size(y)) stop 'Error mixed_central_moment_var_sp: size(x) .ne. size(y)'
1327 3 : if (present(mask)) then
1328 1 : if (size(mask) .ne. size(x)) stop 'Error mixed_central_moment_var_sp: size(mask) .ne. size(x)'
1329 11 : maske = mask
1330 11 : n = real(count(maske), sp)
1331 : else
1332 21 : maske(:) = .true.
1333 2 : n = real(size(x), sp)
1334 : end if
1335 3 : if (n .le. (2.0_sp + tiny(2.0_sp))) stop 'mixed_central_moment_var_sp: n must be at least 3'
1336 :
1337 3 : u2r2s = mixed_central_moment(x, y, 2 * r, 2 * s, mask = maske)
1338 3 : urs = mixed_central_moment(x, y, r, s, mask = maske)
1339 3 : urm1s = mixed_central_moment(x, y, r - 1, s, mask = maske)
1340 3 : u20 = mixed_central_moment(x, y, 2, 0, mask = maske)
1341 3 : urp1s = mixed_central_moment(x, y, r + 1, s, mask = maske)
1342 3 : ursm1 = mixed_central_moment(x, y, r, s - 1, mask = maske)
1343 3 : u02 = mixed_central_moment(x, y, 0, 2, mask = maske)
1344 3 : ursp1 = mixed_central_moment(x, y, r, s + 1, mask = maske)
1345 3 : u11 = mixed_central_moment(x, y, 1, 1, mask = maske)
1346 3 : rr = real(r, sp)
1347 3 : ss = real(s, sp)
1348 :
1349 : mixed_central_moment_var_sp = (u2r2s - urs * urs &
1350 : + rr * rr**u20 * urm1s * urm1s + ss * ss * u02 * ursm1 * ursm1 &
1351 : + 2.0_sp * rr * ss * u11 * urm1s * ursm1 &
1352 3 : - 2.0_sp * rr * urp1s * urm1s - 2.0_sp * ss * ursp1 * ursm1) / n
1353 :
1354 3 : END FUNCTION mixed_central_moment_var_sp
1355 :
1356 : ! ------------------------------------------------------------------
1357 :
1358 48 : SUBROUTINE moment_dp(dat, average, variance, skewness, kurtosis, mean, stddev, absdev, mask, sample)
1359 :
1360 : IMPLICIT NONE
1361 :
1362 : REAL(dp), DIMENSION(:), INTENT(IN) :: dat
1363 : REAL(dp), OPTIONAL, INTENT(OUT) :: average
1364 : REAL(dp), OPTIONAL, INTENT(OUT) :: variance
1365 : REAL(dp), OPTIONAL, INTENT(OUT) :: skewness
1366 : REAL(dp), OPTIONAL, INTENT(OUT) :: kurtosis
1367 : REAL(dp), OPTIONAL, INTENT(OUT) :: mean
1368 : REAL(dp), OPTIONAL, INTENT(OUT) :: stddev
1369 : REAL(dp), OPTIONAL, INTENT(OUT) :: absdev
1370 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
1371 : LOGICAL, OPTIONAL, INTENT(IN) :: sample
1372 :
1373 24 : REAL(dp) :: n, div_n ! divisor depending if sample or population moments
1374 :
1375 9776 : REAL(dp) :: ep, ave, var
1376 19480 : REAL(dp), DIMENSION(size(dat)) :: p, s
1377 24 : LOGICAL, DIMENSION(size(dat)) :: maske
1378 :
1379 24 : if (present(mask)) then
1380 24 : if (size(mask) .ne. size(dat)) stop 'Error moment_dp: size(mask) .ne. size(dat)'
1381 9752 : maske = mask
1382 9752 : n = real(count(maske), sp)
1383 : else
1384 0 : maske(:) = .true.
1385 0 : n = real(size(dat), sp)
1386 : end if
1387 24 : if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'moment_dp: n must be at least 2'
1388 :
1389 : ! set divisor for population (n) or sample (n-1) variance
1390 24 : div_n = n - 1.0_dp
1391 24 : if (present(sample)) then
1392 5 : if (.not. sample) div_n = n
1393 : end if
1394 :
1395 : ! Any optional argument
1396 24 : if (.not. (present(average) .or. present(variance) .or. present(skewness) .or. &
1397 : present(kurtosis) .or. present(mean) .or. present(stddev) .or. present(absdev))) return
1398 : ! Average
1399 9752 : ave = sum(dat(:), mask = maske) / n
1400 24 : if (present(average)) average = ave
1401 24 : if (present(mean)) mean = ave
1402 24 : if (.not. (present(variance) .or. present(skewness) .or. &
1403 : present(kurtosis) .or. present(stddev) .or. present(absdev))) return
1404 : ! Absolute deviation
1405 9752 : s(:) = dat(:) - ave
1406 34 : if (present(absdev)) absdev = sum(abs(s(:)), mask = maske) / n
1407 : ! Variance / Standard deviation
1408 24 : if (.not. (present(variance) .or. present(skewness) .or. &
1409 : present(kurtosis) .or. present(stddev))) return
1410 9752 : ep = sum(s(:), mask = maske)
1411 9752 : p(:) = s(:) * s(:)
1412 9752 : var = sum(p(:), mask = maske)
1413 24 : var = (var - ep * ep / n) / div_n
1414 24 : if (present(variance)) variance = var
1415 : ! Standard deviation
1416 24 : if (present(stddev)) stddev = sqrt(var)
1417 24 : if (.not. (present(skewness) .or. present(kurtosis))) return
1418 : ! Skewness
1419 1 : if (abs(var) .lt. tiny(0.0_dp)) stop 'moment_dp: no skewness or kurtosis when zero variance'
1420 11 : p(:) = p(:) * s(:)
1421 1 : if (present(skewness)) then
1422 11 : skewness = sum(p(:), mask = maske)
1423 1 : skewness = skewness / (n * stddev * stddev * stddev)
1424 : end if
1425 : ! Kurtosis
1426 1 : if (present(kurtosis)) then
1427 11 : p(:) = p(:) * s(:)
1428 11 : kurtosis = sum(p(:), mask = maske)
1429 1 : kurtosis = kurtosis / (n * variance * variance) - 3.0_dp
1430 : end if
1431 :
1432 27 : END SUBROUTINE moment_dp
1433 :
1434 :
1435 48 : SUBROUTINE moment_sp(dat, average, variance, skewness, kurtosis, mean, stddev, absdev, mask, sample)
1436 :
1437 : IMPLICIT NONE
1438 :
1439 : REAL(sp), DIMENSION(:), INTENT(IN) :: dat
1440 : REAL(sp), OPTIONAL, INTENT(OUT) :: average
1441 : REAL(sp), OPTIONAL, INTENT(OUT) :: variance
1442 : REAL(sp), OPTIONAL, INTENT(OUT) :: skewness
1443 : REAL(sp), OPTIONAL, INTENT(OUT) :: kurtosis
1444 : REAL(sp), OPTIONAL, INTENT(OUT) :: mean
1445 : REAL(sp), OPTIONAL, INTENT(OUT) :: stddev
1446 : REAL(sp), OPTIONAL, INTENT(OUT) :: absdev
1447 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
1448 : LOGICAL, OPTIONAL, INTENT(IN) :: sample
1449 :
1450 24 : REAL(sp) :: n, div_n ! divisor depending if sample or population moments
1451 :
1452 9776 : REAL(sp) :: ep, ave, var
1453 19480 : REAL(sp), DIMENSION(size(dat)) :: p, s
1454 24 : LOGICAL, DIMENSION(size(dat)) :: maske
1455 :
1456 24 : if (present(mask)) then
1457 24 : if (size(mask) .ne. size(dat)) stop 'Error moment_sp: size(mask) .ne. size(dat)'
1458 9752 : maske = mask
1459 9752 : n = real(count(maske), sp)
1460 : else
1461 0 : maske(:) = .true.
1462 0 : n = real(size(dat), sp)
1463 : end if
1464 24 : if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'moment_sp: n must be at least 2'
1465 :
1466 : ! set divisor for population (n) or sample (n-1) variance
1467 24 : div_n = n - 1.0_sp
1468 24 : if (present(sample)) then
1469 5 : if (.not. sample) div_n = n
1470 : end if
1471 :
1472 : ! Any optional argument
1473 24 : if (.not. (present(average) .or. present(variance) .or. present(skewness) .or. &
1474 : present(kurtosis) .or. present(mean) .or. present(stddev) .or. present(absdev))) return
1475 : ! Average
1476 9752 : ave = sum(dat(:), mask = maske) / n
1477 24 : if (present(average)) average = ave
1478 24 : if (present(mean)) mean = ave
1479 24 : if (.not. (present(variance) .or. present(skewness) .or. &
1480 : present(kurtosis) .or. present(stddev) .or. present(absdev))) return
1481 : ! Absolute deviation
1482 9752 : s(:) = dat(:) - ave
1483 34 : if (present(absdev)) absdev = sum(abs(s(:)), mask = maske) / n
1484 : ! Variance / Standard deviation
1485 24 : if (.not. (present(variance) .or. present(skewness) .or. &
1486 : present(kurtosis) .or. present(stddev))) return
1487 9752 : ep = sum(s(:), mask = maske)
1488 9752 : p(:) = s(:) * s(:)
1489 9752 : var = sum(p(:), mask = maske)
1490 24 : var = (var - ep * ep / n) / div_n
1491 24 : if (present(variance)) variance = var
1492 : ! Standard deviation
1493 24 : if (present(stddev)) stddev = sqrt(var)
1494 24 : if (.not. (present(skewness) .or. present(kurtosis))) return
1495 : ! Skewness
1496 1 : if (abs(var) .lt. tiny(0.0_sp)) stop 'moment_sp: no skewness or kurtosis when zero variance'
1497 11 : p(:) = p(:) * s(:)
1498 1 : if (present(skewness)) then
1499 11 : skewness = sum(p(:), mask = maske)
1500 1 : skewness = skewness / (n * stddev * stddev * stddev)
1501 : end if
1502 : ! Kurtosis
1503 1 : if (present(kurtosis)) then
1504 11 : p(:) = p(:) * s(:)
1505 11 : kurtosis = sum(p(:), mask = maske)
1506 1 : kurtosis = kurtosis / (n * variance * variance) - 3.0_sp
1507 : end if
1508 :
1509 48 : END SUBROUTINE moment_sp
1510 :
1511 : ! ------------------------------------------------------------------
1512 :
1513 364574 : FUNCTION stddev_dp(dat, mask)
1514 :
1515 : IMPLICIT NONE
1516 :
1517 : REAL(dp), DIMENSION(:), INTENT(IN) :: dat
1518 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
1519 : REAL(dp) :: stddev_dp
1520 :
1521 182287 : REAL(dp) :: n
1522 :
1523 19036075 : REAL(dp) :: ep, ave, var
1524 37525289 : REAL(dp), DIMENSION(size(dat)) :: p, s
1525 182287 : LOGICAL, DIMENSION(size(dat)) :: maske
1526 :
1527 182287 : if (present(mask)) then
1528 33 : if (size(mask) .ne. size(dat)) stop 'Error stddev_dp: size(mask) .ne. size(dat)'
1529 4435 : maske = mask
1530 4435 : n = real(count(maske), dp)
1531 : else
1532 18849353 : maske(:) = .true.
1533 182254 : n = real(size(dat), dp)
1534 : end if
1535 182287 : if (n .le. (1.0_dp + tiny(1.0_dp))) stop 'stddev_dp: n must be at least 2'
1536 :
1537 : ! Average
1538 18853788 : ave = sum(dat(:), mask = maske) / n
1539 18853788 : s(:) = dat(:) - ave
1540 : ! Variance / Standard deviation
1541 18853788 : ep = sum(s(:), mask = maske)
1542 18853788 : p(:) = s(:) * s(:)
1543 18853788 : var = sum(p(:), mask = maske)
1544 182287 : var = (var - ep * ep / n) / (n - 1.0_dp)
1545 182287 : stddev_dp = sqrt(var)
1546 :
1547 24 : END FUNCTION stddev_dp
1548 :
1549 :
1550 70 : FUNCTION stddev_sp(dat, mask)
1551 :
1552 : IMPLICIT NONE
1553 :
1554 : REAL(sp), DIMENSION(:), INTENT(IN) :: dat
1555 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
1556 : REAL(sp) :: stddev_sp
1557 :
1558 35 : REAL(sp) :: n
1559 :
1560 4491 : REAL(sp) :: ep, ave, var
1561 8877 : REAL(sp), DIMENSION(size(dat)) :: p, s
1562 35 : LOGICAL, DIMENSION(size(dat)) :: maske
1563 :
1564 35 : if (present(mask)) then
1565 33 : if (size(mask) .ne. size(dat)) stop 'Error stddev_sp: size(mask) .ne. size(dat)'
1566 4435 : maske = mask
1567 4435 : n = real(count(maske), sp)
1568 : else
1569 21 : maske(:) = .true.
1570 2 : n = real(size(dat), sp)
1571 : end if
1572 35 : if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'stddev_sp: n must be at least 2'
1573 :
1574 : ! Average
1575 4456 : ave = sum(dat(:), mask = maske) / n
1576 4456 : s(:) = dat(:) - ave
1577 : ! Variance / Standard deviation
1578 4456 : ep = sum(s(:), mask = maske)
1579 4456 : p(:) = s(:) * s(:)
1580 4456 : var = sum(p(:), mask = maske)
1581 35 : var = (var - ep * ep / n) / (n - 1.0_sp)
1582 35 : stddev_sp = sqrt(var)
1583 :
1584 182287 : END FUNCTION stddev_sp
1585 :
1586 : ! ------------------------------------------------------------------
1587 :
1588 6 : FUNCTION skewness_dp(dat, mask)
1589 :
1590 : IMPLICIT NONE
1591 :
1592 : REAL(dp), DIMENSION(:), INTENT(IN) :: dat
1593 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
1594 : REAL(dp) :: skewness_dp
1595 :
1596 3 : REAL(dp) :: n
1597 :
1598 35 : REAL(dp) :: ep, ave, var, stddev
1599 61 : REAL(dp), DIMENSION(size(dat)) :: p, s
1600 3 : LOGICAL, DIMENSION(size(dat)) :: maske
1601 :
1602 3 : if (present(mask)) then
1603 1 : if (size(mask) .ne. size(dat)) stop 'Error skewness_dp: size(mask) .ne. size(dat)'
1604 11 : maske = mask
1605 11 : n = real(count(maske), dp)
1606 : else
1607 21 : maske(:) = .true.
1608 2 : n = real(size(dat), dp)
1609 : end if
1610 3 : if (n .le. (1.0_dp + tiny(1.0_dp))) stop 'skewness_dp: n must be at least 2'
1611 :
1612 : ! Average
1613 32 : ave = sum(dat(:), mask = maske) / n
1614 32 : s(:) = dat(:) - ave
1615 : ! Variance / Standard deviation
1616 32 : ep = sum(s(:), mask = maske)
1617 32 : p(:) = s(:) * s(:)
1618 32 : var = sum(p(:), mask = maske)
1619 3 : var = (var - ep * ep / n) / (n - 1.0_dp)
1620 3 : stddev = sqrt(var)
1621 : ! Skewness
1622 3 : if (abs(var) .lt. tiny(0.0_dp)) stop 'skewness_dp: no skewness when zero variance'
1623 32 : p(:) = p(:) * s(:)
1624 32 : skewness_dp = sum(p(:), mask = maske)
1625 3 : skewness_dp = skewness_dp / (n * stddev * stddev * stddev)
1626 :
1627 35 : END FUNCTION skewness_dp
1628 :
1629 :
1630 6 : FUNCTION skewness_sp(dat, mask)
1631 :
1632 : IMPLICIT NONE
1633 :
1634 : REAL(sp), DIMENSION(:), INTENT(IN) :: dat
1635 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
1636 : REAL(sp) :: skewness_sp
1637 :
1638 3 : REAL(sp) :: n
1639 :
1640 35 : REAL(sp) :: ep, ave, var, stddev
1641 61 : REAL(sp), DIMENSION(size(dat)) :: p, s
1642 3 : LOGICAL, DIMENSION(size(dat)) :: maske
1643 :
1644 3 : if (present(mask)) then
1645 1 : if (size(mask) .ne. size(dat)) stop 'Error skewness_sp: size(mask) .ne. size(dat)'
1646 11 : maske = mask
1647 11 : n = real(count(maske), sp)
1648 : else
1649 21 : maske(:) = .true.
1650 2 : n = real(size(dat), sp)
1651 : end if
1652 3 : if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'skewness_sp: n must be at least 2'
1653 :
1654 : ! Average
1655 32 : ave = sum(dat(:), mask = maske) / n
1656 32 : s(:) = dat(:) - ave
1657 : ! Variance / Standard deviation
1658 32 : ep = sum(s(:), mask = maske)
1659 32 : p(:) = s(:) * s(:)
1660 32 : var = sum(p(:), mask = maske)
1661 3 : var = (var - ep * ep / n) / (n - 1.0_sp)
1662 3 : stddev = sqrt(var)
1663 : ! Skewness
1664 3 : if (abs(var) .lt. tiny(0.0_sp)) stop 'skewness_sp: no skewness when zero variance'
1665 32 : p(:) = p(:) * s(:)
1666 32 : skewness_sp = sum(p(:), mask = maske)
1667 3 : skewness_sp = skewness_sp / (n * stddev * stddev * stddev)
1668 :
1669 3 : END FUNCTION skewness_sp
1670 :
1671 : ! ------------------------------------------------------------------
1672 :
1673 6 : FUNCTION variance_dp(dat, mask)
1674 :
1675 : IMPLICIT NONE
1676 :
1677 : REAL(dp), DIMENSION(:), INTENT(IN) :: dat
1678 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
1679 : REAL(dp) :: variance_dp
1680 :
1681 3 : REAL(dp) :: n
1682 :
1683 35 : REAL(dp) :: ep, ave, var
1684 67 : REAL(dp), DIMENSION(size(dat)) :: p, s
1685 3 : LOGICAL, DIMENSION(size(dat)) :: maske
1686 :
1687 3 : if (present(mask)) then
1688 1 : if (size(mask) .ne. size(dat)) stop 'Error variance_dp: size(mask) .ne. size(dat)'
1689 12 : maske = mask
1690 11 : n = real(count(maske), dp)
1691 : else
1692 21 : maske(:) = .true.
1693 2 : n = real(size(dat), dp)
1694 : end if
1695 3 : if (n .le. (1.0_dp + tiny(1.0_dp))) stop 'variance_dp: n must be at least 2'
1696 :
1697 : ! Average
1698 35 : ave = sum(dat(:), mask = maske) / n
1699 32 : s(:) = dat(:) - ave
1700 : ! Variance / Standard deviation
1701 32 : ep = sum(s(:), mask = maske)
1702 32 : p(:) = s(:) * s(:)
1703 32 : var = sum(p(:), mask = maske)
1704 3 : variance_dp = (var - ep * ep / n) / (n - 1.0_dp)
1705 :
1706 3 : END FUNCTION variance_dp
1707 :
1708 :
1709 3 : FUNCTION variance_sp(dat, mask)
1710 :
1711 : IMPLICIT NONE
1712 :
1713 : REAL(sp), DIMENSION(:), INTENT(IN) :: dat
1714 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
1715 : REAL(sp) :: variance_sp
1716 :
1717 3 : REAL(sp) :: n
1718 :
1719 35 : REAL(sp) :: ep, ave, var
1720 67 : REAL(sp), DIMENSION(size(dat)) :: p, s
1721 3 : LOGICAL, DIMENSION(size(dat)) :: maske
1722 :
1723 3 : if (present(mask)) then
1724 1 : if (size(mask) .ne. size(dat)) stop 'Error variance_sp: size(mask) .ne. size(dat)'
1725 12 : maske = mask
1726 11 : n = real(count(maske), sp)
1727 : else
1728 21 : maske(:) = .true.
1729 2 : n = real(size(dat), sp)
1730 : end if
1731 3 : if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'variance_sp: n must be at least 2'
1732 :
1733 : ! Average
1734 35 : ave = sum(dat(:), mask = maske) / n
1735 32 : s(:) = dat(:) - ave
1736 : ! Variance / Standard deviation
1737 32 : ep = sum(s(:), mask = maske)
1738 32 : p(:) = s(:) * s(:)
1739 32 : var = sum(p(:), mask = maske)
1740 3 : variance_sp = (var - ep * ep / n) / (n - 1.0_sp)
1741 :
1742 3 : END FUNCTION variance_sp
1743 :
1744 : ! ------------------------------------------------------------------
1745 :
1746 : END MODULE mo_moment
|
