libsim Versione 7.2.6
simple_stat.f90
1! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2! authors:
3! Davide Cesari <dcesari@arpa.emr.it>
4! Paolo Patruno <ppatruno@arpa.emr.it>
5
6! This program is free software; you can redistribute it and/or
7! modify it under the terms of the GNU General Public License as
8! published by the Free Software Foundation; either version 2 of
9! the License, or (at your option) any later version.
10
11! This program is distributed in the hope that it will be useful,
12! but WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14! GNU General Public License for more details.
15
16! You should have received a copy of the GNU General Public License
17! along with this program. If not, see <http://www.gnu.org/licenses/>.
18!> Module for basic statistical computations taking into account missing data.
19!! The functions provided are average, variance, linear correlation
20!! coefficient and percentile for a population. Missing data are
21!! excluded from computations and from weighing. Both \a REAL and \a
22!! DOUBLE \a PRECISION functions are available.
23!!
24!! \ingroup base
29IMPLICIT NONE
30
31!> Compute the average of the random variable provided,
32!! taking into account missing data. The result is of type \a REAL or
33!! \a DOUBLE \a PRECISION according to the type of \a sample.
34!!
35!! REAL or DOUBLE PRECISION FUNCTION stat_average()
36!! \param sample(:) REAL,INTENT(in) or DOUBLE PRECISION,INTENT(in) the variable to be averaged
37!! \param mask(:) LOGICAL,OPTIONAL,INTENT(in) additional mask to be and'ed with missing values
38!! \param nomiss LOGICAL,OPTIONAL,INTENT(in) if provided and \a .TRUE. it disables all the checks for missing data and empty sample and enables the use of a fast algorithm
39INTERFACE stat_average
40 MODULE PROCEDURE stat_averager, stat_averaged
41END INTERFACE
42
43!> Compute the variance of the random variable provided, taking into
44!! account missing data. The average can be returned as an optional
45!! parameter since it is computed anyway. The result is of type \a REAL or
46!! \a DOUBLE \a PRECISION according to the type of \a sample.
47!!
48!! REAL or DOUBLE PRECISION FUNCTION stat_variance()
49!! \param sample(:) REAL,INTENT(in) or DOUBLE PRECISION,INTENT(in) the variable for which variance has to be computed
50!! \param average REAL,OPTIONAL,INTENT(out) or DOUBLE PRECISION,OPTIONAL,INTENT(out) the average of the variable can optionally be returned
51!! \param mask(:) LOGICAL,OPTIONAL,INTENT(in) additional mask to be and'ed with missing values
52!! \param nomiss LOGICAL,OPTIONAL,INTENT(in) if provided and \a .TRUE. it disables all the checks for missing data and empty sample and enables the use of a fast algorithm
53!! \param nm1 LOGICAL,OPTIONAL,INTENT(in) if provided and \a .TRUE. it computes the variance dividing by \a n - 1 rather than by \a n
55 MODULE PROCEDURE stat_variancer, stat_varianced
56END INTERFACE
57
58!> Compute the standard deviation of the random variable provided, taking into
59!! account missing data. The average can be returned as an optional
60!! parameter since it is computed anyway. The result is of type \a REAL or
61!! \a DOUBLE \a PRECISION according to the type of \a sample.
62!!
63!! REAL or DOUBLE PRECISION FUNCTION stat_stddev()
64!! \param sample(:) REAL,INTENT(in) or DOUBLE PRECISION,INTENT(in) the variable for which standard deviation has to be computed
65!! \param average REAL,OPTIONAL,INTENT(out) or DOUBLE PRECISION,OPTIONAL,INTENT(out) the average of the variable can optionally be returned
66!! \param mask(:) LOGICAL,OPTIONAL,INTENT(in) additional mask to be and'ed with missing values
67!! \param nomiss LOGICAL,OPTIONAL,INTENT(in) if provided and \a .TRUE. it disables all the checks for missing data and empty sample and enables the use of a fast algorithm
68!! \param nm1 LOGICAL,OPTIONAL,INTENT(in) if provided and \a .TRUE. it computes the variance dividing by \a n - 1 rather than by \a n
69INTERFACE stat_stddev
70 MODULE PROCEDURE stat_stddevr, stat_stddevd
71END INTERFACE
72
73!> Compute the linear correlation coefficient between the two random
74!! variables provided, taking into account missing data. Data are
75!! considered missing when at least one variable has a missing value.
76!! The average and the variance of each variable can be returned as
77!! optional parameters since they are computed anyway. The result is
78!! of type \a REAL or \a DOUBLE \a PRECISION according to the type of
79!! \a sample1 and \a sample2.
80!!
81!! REAL or DOUBLE PRECISION FUNCTION stat_linear_corr()
82!! \param sample1(:) REAL,INTENT(in) or DOUBLE PRECISION,INTENT(in) the first variable
83!! \param sample2(:) REAL,INTENT(in) or DOUBLE PRECISION,INTENT(in) the second variable
84!! \param average1 REAL,OPTIONAL,INTENT(out) or DOUBLE PRECISION,OPTIONAL,INTENT(out) the average of the first variable can optionally be returned
85!! \param average2 REAL,OPTIONAL,INTENT(out) or DOUBLE PRECISION,OPTIONAL,INTENT(out) the average of the second variable can optionally be returned
86!! \param variance1 REAL,OPTIONAL,INTENT(out) or DOUBLE PRECISION,OPTIONAL,INTENT(out) the variance of the first variable can optionally be returned
87!! \param variance2 REAL,OPTIONAL,INTENT(out) or DOUBLE PRECISION,OPTIONAL,INTENT(out) the variance of the second variable can optionally be returned
88!! \param mask(:) LOGICAL,OPTIONAL,INTENT(in) additional mask to be and'ed with missing values
90 MODULE PROCEDURE stat_linear_corrr, stat_linear_corrd
91END INTERFACE
92
93!> Compute the linear regression coefficients between the two random
94!! variables provided, taking into account missing data. Data are
95!! considered missing when at least one variable has a missing value.
96!! The regression is computed using the method of linear least
97!! squares. The input and output parameters are either \a REAL or \a
98!! DOUBLE \a PRECISION.
99!!
100!! SUBROUTINE stat_linear_regression()
101!! \param sample1(:) REAL,INTENT(in) or DOUBLE PRECISION,INTENT(in) the first variable
102!! \param sample2(:) REAL,INTENT(in) or DOUBLE PRECISION,INTENT(in) the second variable
103!! \param alpha0 REAL,INTENT(out) or DOUBLE PRECISION,INTENT(out) the 0-th order coefficient of the regression computed
104!! \param alpha1 REAL,INTENT(out) or DOUBLE PRECISION,INTENT(out) the first order coefficient of the regression computed
105!! \param mask(:) LOGICAL,OPTIONAL,INTENT(in) additional mask to be and'ed with missing values
106!! \param nomiss LOGICAL,OPTIONAL,INTENT(in) if provided and \a .TRUE. it disables all the checks for missing data and empty sample and enables the use of a fast algorithm
108 MODULE PROCEDURE stat_linear_regressionr, stat_linear_regressiond
109END INTERFACE
110
111!> Compute a set of percentiles for a random variable.
112!! The result is a 1-d array of the same size as the array of
113!! percentile values provided. Percentiles are computed by linear
114!! interpolation. The result is of type \a REAL or \a DOUBLE \a
115!! PRECISION according to the type of \a sample.
116!!
117!! REAL or DOUBLE PRECISION FUNCTION stat_percentile()
118!! \param sample(:) REAL,INTENT(in) or DOUBLE PRECISION,INTENT(in) the variable for which percentiles have to be computed
119!! \param perc_vals(:) REAL,INTENT(in) or DOUBLE PRECISION,INTENT(in) the percentiles values to be computed, between 0. and 100.
120!! \param mask(:) LOGICAL,OPTIONAL,INTENT(in) additional mask to be and'ed with missing values
122 MODULE PROCEDURE stat_percentiler, stat_percentiled
123END INTERFACE
124
125!> Bin a sample into equally spaced intervals to form a histogram.
126!! The sample is binned into the requested number of intervals (bins)
127!! and the population of each bin is returned in the allocatable array
128!! \a bin(:). If the operation cannot be performed \a bin is not
129!! allocated. If the lower and upper bound of the histogram are not
130!! provided, they are computed as the minimum and maximum of the
131!! sample.
132!!
133!! SUBROUTINE stat_bin()
134!! \param sample(:) REAL or DOUBLE PRECISION,INTENT(in) the variable to be binned
135!! \param bin(:) INTEGER,INTENT(out),ALLOCATABLE the array with the population of each bin, dimensioned as (nbin)
136!! \param nbin INTEGER,INTENT(in) the number of bins requested
137!! \param start REAL or DOUBLE PRECISION,INTENT(in),OPTIONAL the start of the overall histogram interval
138!! \param finish REAL or DOUBLE PRECISION,INTENT(in),OPTIONAL the end of the overall histogram interval
139!! \param mask(:) LOGICAL,OPTIONAL,INTENT(in) additional mask to be and'ed with missing values
140!! \param binbounds(:) REAL or DOUBLE PRECISION,INTENT(out),ALLOCATABLE,OPTIONAL the boundary of each bin, dimensioned as (nbin+1)
141INTERFACE stat_bin
142 MODULE PROCEDURE stat_binr, stat_bind
143END INTERFACE stat_bin
144
145!> Compute the mode of the random variable provided
146!! taking into account missing data.
147!! The mode is computed by binning the sample into the requested
148!! number of intervals (bins) and the mode is computed as the midpoint
149!! of the most populated bin; in case of multi-modality, the smallest
150!! value is returned. If the operation cannot be performed a missing
151!! value is returned. See the description of simple_stat::stat_bin
152!! subroutine for details about the binning procedure.
153!!
154!! REAL or DOUBLE PRECISION FUNCTION stat_mode_histogram()
155!! \param sample(:) REAL or DOUBLE PRECISION,INTENT(in) the variable for which the mode has to be computed
156!! \param nbin INTEGER,INTENT(in) the number of bins requested
157!! \param start REAL or DOUBLE PRECISION,INTENT(in),OPTIONAL the start of the overall histogram interval
158!! \param finish REAL or DOUBLE PRECISION,INTENT(in),OPTIONAL the end of the overall histogram interval
159!! \param mask(:) LOGICAL,OPTIONAL,INTENT(in) additional mask to be and'ed with missing values
161 MODULE PROCEDURE stat_mode_histogramr, stat_mode_histogramd
162END INTERFACE stat_mode_histogram
163
164PRIVATE
167 normalizeddensityindex
168
169CONTAINS
170
171
172FUNCTION stat_averager(sample, mask, nomiss) RESULT(average)
173REAL,INTENT(in) :: sample(:)
174LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
175LOGICAL,OPTIONAL,INTENT(in) :: nomiss ! informs that the sample does not contain missing data and is not of zero length (for speedup)
176
177REAL :: average
178
179INTEGER :: sample_count
180LOGICAL :: sample_mask(SIZE(sample))
181
182IF (optio_log(nomiss)) THEN
183 average = sum(sample)/SIZE(sample)
184ELSE
185 sample_mask = (sample /= rmiss)
186 IF (PRESENT(mask)) sample_mask = sample_mask .AND. mask
187 sample_count = count(sample_mask)
188
189 IF (sample_count > 0) THEN
190! compute average
191 average = sum(sample, mask=sample_mask)/sample_count
192 ELSE
193 average = rmiss
194 ENDIF
195ENDIF
196
197END FUNCTION stat_averager
198
199
200FUNCTION stat_averaged(sample, mask, nomiss) RESULT(average)
201DOUBLE PRECISION,INTENT(in) :: sample(:)
202LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
203LOGICAL,OPTIONAL,INTENT(in) :: nomiss ! informs that the sample does not contain missing data and is not of zero length (for speedup)
204
205DOUBLE PRECISION :: average
206
207INTEGER :: sample_count
208LOGICAL :: sample_mask(SIZE(sample))
209
210IF (optio_log(nomiss)) THEN
211 average = sum(sample)/SIZE(sample)
212ELSE
213 sample_mask = (sample /= dmiss)
214 IF (PRESENT(mask)) sample_mask = sample_mask .AND. mask
215 sample_count = count(sample_mask)
216
217 IF (sample_count > 0) THEN
218! compute average
219 average = sum(sample, mask=sample_mask)/sample_count
220 ELSE
221 average = dmiss
222 ENDIF
223ENDIF
224
225END FUNCTION stat_averaged
226
227
228FUNCTION stat_variancer(sample, average, mask, nomiss, nm1) RESULT(variance)
229REAL,INTENT(in) :: sample(:)
230REAL,OPTIONAL,INTENT(out) :: average
231LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
232LOGICAL,OPTIONAL,INTENT(in) :: nomiss
233LOGICAL,OPTIONAL,INTENT(in) :: nm1
234
235REAL :: variance
236
237REAL :: laverage
238INTEGER :: sample_count, i
239LOGICAL :: sample_mask(SIZE(sample))
240
241IF (optio_log(nomiss)) THEN
242! compute average
243 laverage = sum(sample)/SIZE(sample)
244 IF (PRESENT(average)) average = laverage
245 IF (optio_log(nm1)) THEN
246 variance = sum((sample-laverage)**2)/max(SIZE(sample)-1,1)
247 ELSE
248 variance = sum((sample-laverage)**2)/SIZE(sample)
249 ENDIF
250
251ELSE
252 sample_mask = (sample /= rmiss)
253 IF (PRESENT(mask)) sample_mask = sample_mask .AND. mask
254 sample_count = count(sample_mask)
255
256 IF (sample_count > 0) THEN
257! compute average
258 laverage = sum(sample, mask=sample_mask)/sample_count
259 IF (PRESENT(average)) average = laverage
260! compute sum of squares and variance
261 variance = 0.
262 DO i = 1, SIZE(sample)
263 IF (sample_mask(i)) variance = variance + (sample(i)-laverage)**2
264 ENDDO
265 IF (optio_log(nm1)) THEN
266 variance = variance/max(sample_count-1,1)
267 ELSE
268 variance = variance/sample_count
269 ENDIF
270 ELSE
271 IF (PRESENT(average)) average = rmiss
272 variance = rmiss
273 ENDIF
274
275ENDIF
276
277END FUNCTION stat_variancer
278
279
280FUNCTION stat_varianced(sample, average, mask, nomiss, nm1) RESULT(variance)
281DOUBLE PRECISION,INTENT(in) :: sample(:)
282DOUBLE PRECISION,OPTIONAL,INTENT(out) :: average
283LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
284LOGICAL,OPTIONAL,INTENT(in) :: nomiss
285LOGICAL,OPTIONAL,INTENT(in) :: nm1
286
287DOUBLE PRECISION :: variance
288
289DOUBLE PRECISION :: laverage
290INTEGER :: sample_count, i
291LOGICAL :: sample_mask(SIZE(sample))
292
293IF (optio_log(nomiss)) THEN
294! compute average
295 laverage = sum(sample)/SIZE(sample)
296 IF (PRESENT(average)) average = laverage
297 IF (optio_log(nm1)) THEN
298 variance = sum((sample-laverage)**2)/max(SIZE(sample)-1,1)
299 ELSE
300 variance = sum((sample-laverage)**2)/SIZE(sample)
301 ENDIF
302
303ELSE
304 sample_mask = (sample /= dmiss)
305 IF (PRESENT(mask)) sample_mask = sample_mask .AND. mask
306 sample_count = count(sample_mask)
307
308 IF (sample_count > 0) THEN
309! compute average
310 laverage = sum(sample, mask=sample_mask)/sample_count
311 IF (PRESENT(average)) average = laverage
312! compute sum of squares and variance
313 variance = 0.
314 DO i = 1, SIZE(sample)
315 IF (sample_mask(i)) variance = variance + (sample(i)-laverage)**2
316 ENDDO
317 IF (optio_log(nm1)) THEN
318 variance = variance/max(sample_count-1,1)
319 ELSE
320 variance = variance/sample_count
321 ENDIF
322 ELSE
323 IF (PRESENT(average)) average = dmiss
324 variance = dmiss
325 ENDIF
326
327ENDIF
328
329END FUNCTION stat_varianced
330
331
332FUNCTION stat_stddevr(sample, average, mask, nomiss, nm1) RESULT(stddev)
333REAL,INTENT(in) :: sample(:)
334REAL,OPTIONAL,INTENT(out) :: average
335LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
336LOGICAL,OPTIONAL,INTENT(in) :: nomiss
337LOGICAL,OPTIONAL,INTENT(in) :: nm1
338
339REAL :: stddev
340
341stddev = stat_variance(sample, average, mask, nomiss, nm1)
342IF (c_e(stddev)) stddev = sqrt(stddev)
343
344END FUNCTION stat_stddevr
345
346
347FUNCTION stat_stddevd(sample, average, mask, nomiss, nm1) RESULT(stddev)
348DOUBLE PRECISION,INTENT(in) :: sample(:)
349DOUBLE PRECISION,OPTIONAL,INTENT(out) :: average
350LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
351LOGICAL,OPTIONAL,INTENT(in) :: nomiss
352LOGICAL,OPTIONAL,INTENT(in) :: nm1
353
354DOUBLE PRECISION :: stddev
355
356stddev = stat_variance(sample, average, mask, nomiss, nm1)
357IF (c_e(stddev)) stddev = sqrt(stddev)
358
359END FUNCTION stat_stddevd
360
361
362FUNCTION stat_linear_corrr(sample1, sample2, average1, average2, &
363 variance1, variance2, mask, nomiss) RESULT(linear_corr)
364REAL,INTENT(in) :: sample1(:)
365REAL,INTENT(in) :: sample2(:)
366REAL,OPTIONAL,INTENT(out) :: average1
367REAL,OPTIONAL,INTENT(out) :: average2
368REAL,OPTIONAL,INTENT(out) :: variance1
369REAL,OPTIONAL,INTENT(out) :: variance2
370LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
371LOGICAL,OPTIONAL,INTENT(in) :: nomiss
372
373REAL :: linear_corr
374
375REAL :: laverage1, laverage2, lvariance1, lvariance2
376INTEGER :: sample_count, i
377LOGICAL :: sample_mask(SIZE(sample1))
378
379IF (SIZE(sample1) /= SIZE(sample2)) THEN
380 IF (PRESENT(average1)) average1 = rmiss
381 IF (PRESENT(average2)) average2 = rmiss
382 IF (PRESENT(variance1)) variance1 = rmiss
383 IF (PRESENT(variance2)) variance2 = rmiss
384 linear_corr = rmiss
385 RETURN
386ENDIF
387
388sample_mask = (sample1 /= rmiss .AND. sample2 /= rmiss)
389IF (PRESENT(mask)) sample_mask = sample_mask .AND. mask
390sample_count = count(sample_mask)
391IF (sample_count > 0) THEN
392! compute averages
393 laverage1 = sum(sample1, mask=sample_mask)/sample_count
394 laverage2 = sum(sample2, mask=sample_mask)/sample_count
395 IF (PRESENT(average1)) average1 = laverage1
396 IF (PRESENT(average2)) average2 = laverage2
397! compute sum of squares and variances
398 lvariance1 = 0.
399 lvariance2 = 0.
400 DO i = 1, SIZE(sample1)
401 IF (sample_mask(i))THEN
402 lvariance1 = lvariance1 + (sample1(i)-laverage1)**2
403 lvariance2 = lvariance2 + (sample2(i)-laverage2)**2
404 ENDIF
405 ENDDO
406 lvariance1 = lvariance1/sample_count
407 lvariance2 = lvariance2/sample_count
408 IF (PRESENT(variance1)) variance1 = lvariance1
409 IF (PRESENT(variance2)) variance2 = lvariance2
410! compute correlation
411 linear_corr = 0.
412 DO i = 1, SIZE(sample1)
413 IF (sample_mask(i)) linear_corr = linear_corr + &
414 (sample1(i)-laverage1)*(sample2(i)-laverage2)
415 ENDDO
416 linear_corr = linear_corr/sample_count / sqrt(lvariance1*lvariance2)
417ELSE
418 IF (PRESENT(average1)) average1 = rmiss
419 IF (PRESENT(average2)) average2 = rmiss
420 IF (PRESENT(variance1)) variance1 = rmiss
421 IF (PRESENT(variance2)) variance2 = rmiss
422 linear_corr = rmiss
423ENDIF
424
425END FUNCTION stat_linear_corrr
426
427
428FUNCTION stat_linear_corrd(sample1, sample2, average1, average2, &
429 variance1, variance2, mask, nomiss) RESULT(linear_corr)
430DOUBLE PRECISION, INTENT(in) :: sample1(:)
431DOUBLE PRECISION, INTENT(in) :: sample2(:)
432DOUBLE PRECISION, OPTIONAL, INTENT(out) :: average1
433DOUBLE PRECISION, OPTIONAL, INTENT(out) :: average2
434DOUBLE PRECISION, OPTIONAL, INTENT(out) :: variance1
435DOUBLE PRECISION, OPTIONAL, INTENT(out) :: variance2
436LOGICAL, OPTIONAL, INTENT(in) :: mask(:)
437LOGICAL, OPTIONAL, INTENT(in) :: nomiss
438
439DOUBLE PRECISION :: linear_corr
440
441DOUBLE PRECISION :: laverage1, laverage2, lvariance1, lvariance2
442INTEGER :: sample_count, i
443LOGICAL :: sample_mask(SIZE(sample1))
444
445IF (SIZE(sample1) /= SIZE(sample2)) THEN
446 IF (PRESENT(average1)) average1 = dmiss
447 IF (PRESENT(average2)) average2 = dmiss
448 IF (PRESENT(variance1)) variance1 = dmiss
449 IF (PRESENT(variance2)) variance2 = dmiss
450 linear_corr = dmiss
451 RETURN
452ENDIF
453
454sample_mask = (sample1 /= dmiss .AND. sample2 /= dmiss)
455IF (PRESENT(mask)) sample_mask = sample_mask .AND. mask
456sample_count = count(sample_mask)
457IF (sample_count > 0) THEN
458! compute averages
459 laverage1 = sum(sample1, mask=sample_mask)/sample_count
460 laverage2 = sum(sample2, mask=sample_mask)/sample_count
461 IF (PRESENT(average1)) average1 = laverage1
462 IF (PRESENT(average2)) average2 = laverage2
463! compute sum of squares and variances
464 lvariance1 = 0.
465 lvariance2 = 0.
466 DO i = 1, SIZE(sample1)
467 IF (sample_mask(i))THEN
468 lvariance1 = lvariance1 + (sample1(i)-laverage1)**2
469 lvariance2 = lvariance2 + (sample2(i)-laverage2)**2
470 ENDIF
471 ENDDO
472 lvariance1 = lvariance1/sample_count
473 lvariance2 = lvariance2/sample_count
474 IF (PRESENT(variance1)) variance1 = lvariance1
475 IF (PRESENT(variance2)) variance2 = lvariance2
476! compute correlation
477 linear_corr = 0.0d0
478 DO i = 1, SIZE(sample1)
479 IF (sample_mask(i)) linear_corr = linear_corr + &
480 (sample1(i)-laverage1)*(sample2(i)-laverage2)
481 ENDDO
482 linear_corr = linear_corr/sample_count / sqrt(lvariance1*lvariance2)
483ELSE
484 IF (PRESENT(average1)) average1 = dmiss
485 IF (PRESENT(average2)) average2 = dmiss
486 IF (PRESENT(variance1)) variance1 = dmiss
487 IF (PRESENT(variance2)) variance2 = dmiss
488 linear_corr = dmiss
489ENDIF
490
491END FUNCTION stat_linear_corrd
492
493
494SUBROUTINE stat_linear_regressionr(sample1, sample2, alpha0, alpha1, mask)
495REAL,INTENT(in) :: sample1(:)
496REAL,INTENT(in) :: sample2(:)
497REAL,INTENT(out) :: alpha0
498REAL,INTENT(out) :: alpha1
499LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
500
501REAL :: laverage1, laverage2
502INTEGER :: sample_count
503LOGICAL :: sample_mask(SIZE(sample1))
504
505IF (SIZE(sample1) /= SIZE(sample2)) THEN
506 alpha0 = rmiss
507 alpha1 = rmiss
508 RETURN
509ENDIF
510
511sample_mask = (sample1 /= rmiss .AND. sample2 /= rmiss)
512IF (PRESENT(mask)) sample_mask = sample_mask .AND. mask
513sample_count = count(sample_mask)
514
515IF (sample_count > 0) THEN
516 laverage1 = sum(sample1, mask=sample_mask)/sample_count
517 laverage2 = sum(sample2, mask=sample_mask)/sample_count
518 alpha1 = sum((sample1-laverage1)*(sample2-laverage2), mask=sample_mask)/ &
519 sum((sample1-laverage1)**2, mask=sample_mask)
520 alpha0 = laverage1 - alpha1*laverage2
521
522ELSE
523 alpha0 = rmiss
524 alpha1 = rmiss
525ENDIF
526
527END SUBROUTINE stat_linear_regressionr
528
529
530SUBROUTINE stat_linear_regressiond(sample1, sample2, alpha0, alpha1, mask)
531DOUBLE PRECISION,INTENT(in) :: sample1(:)
532DOUBLE PRECISION,INTENT(in) :: sample2(:)
533DOUBLE PRECISION,INTENT(out) :: alpha0
534DOUBLE PRECISION,INTENT(out) :: alpha1
535LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
536
537DOUBLE PRECISION :: laverage1, laverage2
538INTEGER :: sample_count
539LOGICAL :: sample_mask(SIZE(sample1))
540
541IF (SIZE(sample1) /= SIZE(sample2)) THEN
542 alpha0 = dmiss
543 alpha1 = dmiss
544 RETURN
545ENDIF
546
547sample_mask = (sample1 /= dmiss .AND. sample2 /= dmiss)
548IF (PRESENT(mask)) sample_mask = sample_mask .AND. mask
549sample_count = count(sample_mask)
550
551IF (sample_count > 0) THEN
552 laverage1 = sum(sample1, mask=sample_mask)/sample_count
553 laverage2 = sum(sample2, mask=sample_mask)/sample_count
554 alpha1 = sum((sample1-laverage1)*(sample2-laverage2), mask=sample_mask)/ &
555 sum((sample1-laverage1)**2, mask=sample_mask)
556 alpha0 = laverage1 - alpha1*laverage2
557
558ELSE
559 alpha0 = dmiss
560 alpha1 = dmiss
561ENDIF
562
563END SUBROUTINE stat_linear_regressiond
564
565
566FUNCTION stat_percentiler(sample, perc_vals, mask, nomiss) RESULT(percentile)
567REAL,INTENT(in) :: sample(:)
568REAL,INTENT(in) :: perc_vals(:)
569LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
570LOGICAL,OPTIONAL,INTENT(in) :: nomiss
571REAL :: percentile(SIZE(perc_vals))
572
573REAL :: lsample(SIZE(sample)), rindex
574INTEGER :: sample_count, j, iindex
575LOGICAL :: sample_mask(SIZE(sample))
576
577percentile(:) = rmiss
578IF (.NOT.optio_log(nomiss)) THEN
579 sample_mask = c_e(sample)
580 IF (PRESENT(mask)) sample_mask = sample_mask .AND. mask
581 sample_count = count(sample_mask)
582 IF (sample_count == 0) RETURN ! special case
583ELSE
584 sample_count = SIZE(sample)
585ENDIF
586
587IF (sample_count == SIZE(sample)) THEN
588 lsample = sample
589ELSE
590 lsample(1:sample_count) = pack(sample, mask=sample_mask)
591ENDIF
592
593IF (sample_count == 1) THEN ! other special case
594 percentile(:) = lsample(1)
595 RETURN
596ENDIF
597
598
599! this sort is very fast but with a lot of equal values it is very slow and fails
600CALL sort(lsample(1:sample_count))
601
602! this sort is very very slow
603!!$DO j = 2, sample_count
604!!$ v = lsample(j)
605!!$ DO i = j-1, 1, -1
606!!$ IF (v >= lsample(i)) EXIT
607!!$ lsample(i+1) = lsample(i)
608!!$ ENDDO
609!!$ lsample(i+1) = v
610!!$ENDDO
611
612DO j = 1, SIZE(perc_vals)
613 IF (perc_vals(j) >= 0. .AND. perc_vals(j) <= 100.) THEN
614! compute real index of requested percentile in sample
615 rindex = real(sample_count-1, kind=kind(rindex))*perc_vals(j)/100.+1.
616! compute integer index of previous element in sample, beware of corner cases
617 iindex = min(max(int(rindex), 1), sample_count-1)
618! compute linearly interpolated percentile
619
620 percentile(j) = lsample(iindex)*(real(iindex+1, kind=kind(rindex))-rindex) &
621 + lsample(iindex+1)*(rindex-real(iindex, kind=kind(rindex)))
622
623 ENDIF
624ENDDO
625
626END FUNCTION stat_percentiler
627
628
629FUNCTION stat_percentiled(sample, perc_vals, mask, nomiss) RESULT(percentile)
630DOUBLE PRECISION, INTENT(in) :: sample(:)
631DOUBLE PRECISION, INTENT(in) :: perc_vals(:)
632LOGICAL, OPTIONAL, INTENT(in) :: mask(:)
633LOGICAL, OPTIONAL, INTENT(in) :: nomiss
634DOUBLE PRECISION :: percentile(SIZE(perc_vals))
635
636DOUBLE PRECISION :: lsample(SIZE(sample)), rindex
637INTEGER :: sample_count, j, iindex
638LOGICAL :: sample_mask(SIZE(sample))
639
640
641percentile(:) = dmiss
642IF (.NOT.optio_log(nomiss)) THEN
643 sample_mask = (sample /= dmiss)
644 IF (PRESENT(mask)) sample_mask = sample_mask .AND. mask
645 sample_count = count(sample_mask)
646 IF (sample_count == 0) RETURN ! particular case
647ELSE
648 sample_count = SIZE(sample)
649ENDIF
650
651IF (sample_count == SIZE(sample)) THEN
652 lsample = sample
653ELSE
654 lsample(1:sample_count) = pack(sample, mask=sample_mask)
655ENDIF
656
657IF (sample_count == 1) THEN ! other particular case
658 percentile(:) = lsample(1)
659 RETURN
660ENDIF
661
662! this sort is very fast but with a lot of equal values it is very slow and fails
663CALL sort(lsample(1:sample_count))
664
665! this sort is very very slow
666!DO j = 2, sample_count
667! v = lsample(j)
668! DO i = j-1, 1, -1
669! IF (v >= lsample(i)) EXIT
670! lsample(i+1) = lsample(i)
671! ENDDO
672! lsample(i+1) = v
673!ENDDO
674
675DO j = 1, SIZE(perc_vals)
676 IF (perc_vals(j) >= 0.d0 .AND. perc_vals(j) <= 100.d0) THEN
677! compute real index of requested percentile in sample
678 rindex = real(sample_count-1, kind=kind(rindex))*perc_vals(j)/100.d0+1.d0
679! compute integer index of previous element in sample, beware of corner cases
680 iindex = min(max(int(rindex), 1), sample_count-1)
681! compute linearly interpolated percentile
682 percentile(j) = lsample(iindex)*(real(iindex+1, kind=kind(rindex))-rindex) &
683 + lsample(iindex+1)*(rindex-real(iindex, kind=kind(rindex)))
684 ENDIF
685ENDDO
686
687END FUNCTION stat_percentiled
688
689
690SUBROUTINE stat_binr(sample, bin, nbin, start, finish, mask, binbounds)
691REAL,INTENT(in) :: sample(:)
692INTEGER,INTENT(out),ALLOCATABLE :: bin(:)
693INTEGER,INTENT(in) :: nbin
694REAL,INTENT(in),OPTIONAL :: start
695REAL,INTENT(in),OPTIONAL :: finish
696LOGICAL,INTENT(in),OPTIONAL :: mask(:)
697REAL,INTENT(out),ALLOCATABLE,OPTIONAL :: binbounds(:)
698
699INTEGER :: i, ind
700REAL :: lstart, lfinish, incr
701REAL,ALLOCATABLE :: lbinbounds(:)
702LOGICAL :: lmask(SIZE(sample))
703
704! safety checks
705IF (nbin < 1) RETURN
706IF (PRESENT(mask)) THEN
707 lmask = mask
708ELSE
709 lmask =.true.
710ENDIF
711lmask = lmask .AND. c_e(sample)
712IF (count(lmask) < 1) RETURN
713
714lstart = optio_r(start)
715IF (.NOT.c_e(lstart)) lstart = minval(sample, mask=lmask)
716lfinish = optio_r(finish)
717IF (.NOT.c_e(lfinish)) lfinish = maxval(sample, mask=lmask)
718IF (lfinish <= lstart) RETURN
719
720incr = (lfinish-lstart)/nbin
721ALLOCATE(bin(nbin))
722
723ALLOCATE(lbinbounds(nbin+1))
724
725DO i = 1, nbin
726 lbinbounds(i) = lstart + (i-1)*incr
727ENDDO
728lbinbounds(nbin+1) = lfinish ! set exact value to avoid truncation error
729
730DO i = 1, nbin-1
731 bin(i) = count(sample >= lbinbounds(i) .AND. sample < lbinbounds(i+1) .AND. lmask)
732!, .AND. c_e(sample))
733ENDDO
734bin(nbin) = count(sample >= lbinbounds(nbin) .AND. sample <= lbinbounds(nbin+1) .AND. lmask)
735!, .AND. c_e(sample)) ! special case, include upper limit
736
737IF (PRESENT(binbounds)) binbounds = lbinbounds
738
739END SUBROUTINE stat_binr
740
741
742! alternative algorithm, probably faster with big samples, to be tested
743SUBROUTINE stat_binr2(sample, bin, nbin, start, finish, mask, binbounds)
744REAL,INTENT(in) :: sample(:)
745INTEGER,INTENT(out),ALLOCATABLE :: bin(:)
746INTEGER,INTENT(in) :: nbin
747REAL,INTENT(in),OPTIONAL :: start
748REAL,INTENT(in),OPTIONAL :: finish
749LOGICAL,INTENT(in),OPTIONAL :: mask(:)
750REAL,INTENT(out),ALLOCATABLE,OPTIONAL :: binbounds(:)
751
752INTEGER :: i, ind
753REAL :: lstart, lfinish, incr
754LOGICAL :: lmask(SIZE(sample))
755
756! safety checks
757IF (nbin < 1) RETURN
758IF (PRESENT(mask)) THEN
759 lmask = mask
760ELSE
761 lmask =.true.
762ENDIF
763lmask = lmask .AND. c_e(sample)
764IF (count(lmask) < 1) RETURN
765
766lstart = optio_r(start)
767IF (.NOT.c_e(lstart)) lstart = minval(sample, mask=c_e(sample))
768lfinish = optio_r(finish)
769IF (.NOT.c_e(lfinish)) lfinish = maxval(sample, mask=c_e(sample))
770IF (lfinish <= lstart) RETURN
771
772incr = (lfinish-lstart)/nbin
773ALLOCATE(bin(nbin))
774
775bin(:) = 0
776
777DO i = 1, SIZE(sample)
778 IF (lmask(i)) THEN
779 ind = int((sample(i)-lstart)/incr) + 1
780 IF (ind > 0 .AND. ind <= nbin) THEN
781 bin(ind) = bin(ind) + 1
782 ELSE
783 IF (sample(i) == finish) bin(nbin) = bin(nbin) + 1 ! special case, include upper limit
784 ENDIF
785 ENDIF
786ENDDO
787
788IF (PRESENT(binbounds)) THEN
789 ALLOCATE(binbounds(nbin+1))
790 DO i = 1, nbin
791 binbounds(i) = lstart + (i-1)*incr
792 ENDDO
793 binbounds(nbin+1) = lfinish ! set exact value to avoid truncation error
794ENDIF
795
796END SUBROUTINE stat_binr2
797
798
799SUBROUTINE stat_bind(sample, bin, nbin, start, finish, mask, binbounds)
800DOUBLE PRECISION,INTENT(in) :: sample(:)
801INTEGER,INTENT(out),ALLOCATABLE :: bin(:)
802INTEGER,INTENT(in) :: nbin
803DOUBLE PRECISION,INTENT(in),OPTIONAL :: start
804DOUBLE PRECISION,INTENT(in),OPTIONAL :: finish
805LOGICAL,INTENT(in),OPTIONAL :: mask(:)
806DOUBLE PRECISION,INTENT(out),ALLOCATABLE,OPTIONAL :: binbounds(:)
807
808INTEGER :: i, ind
809DOUBLE PRECISION :: lstart, lfinish, incr
810DOUBLE PRECISION,ALLOCATABLE :: lbinbounds(:)
811LOGICAL :: lmask(SIZE(sample))
812
813! safety checks
814IF (nbin < 1) RETURN
815IF (PRESENT(mask)) THEN
816 lmask = mask
817ELSE
818 lmask =.true.
819ENDIF
820lmask = lmask .AND. c_e(sample)
821IF (count(lmask) < 1) RETURN
822
823lstart = optio_d(start)
824IF (.NOT.c_e(lstart)) lstart = minval(sample, mask=lmask)
825lfinish = optio_d(finish)
826IF (.NOT.c_e(lfinish)) lfinish = maxval(sample, mask=lmask)
827IF (lfinish <= lstart) RETURN
828
829incr = (lfinish-lstart)/nbin
830ALLOCATE(bin(nbin))
831
832ALLOCATE(lbinbounds(nbin+1))
833
834DO i = 1, nbin
835 lbinbounds(i) = lstart + (i-1)*incr
836ENDDO
837lbinbounds(nbin+1) = lfinish ! set exact value to avoid truncation error
838
839DO i = 1, nbin-1
840 bin(i) = count(sample >= lbinbounds(i) .AND. sample < lbinbounds(i+1) .AND. lmask)
841!, .AND. c_e(sample))
842ENDDO
843bin(nbin) = count(sample >= lbinbounds(nbin) .AND. sample <= lbinbounds(nbin+1) .AND. lmask)
844!, .AND. c_e(sample)) ! special case, include upper limit
845
846IF (PRESENT(binbounds)) binbounds = lbinbounds
847
848END SUBROUTINE stat_bind
849
850
851FUNCTION stat_mode_histogramr(sample, nbin, start, finish, mask) RESULT(mode)
852REAL,INTENT(in) :: sample(:)
853INTEGER,INTENT(in) :: nbin
854REAL,INTENT(in),OPTIONAL :: start
855REAL,INTENT(in),OPTIONAL :: finish
856LOGICAL,INTENT(in),OPTIONAL :: mask(:)
857
858REAL :: mode
859
860INTEGER :: loc(1)
861INTEGER,ALLOCATABLE :: bin(:)
862REAL,ALLOCATABLE :: binbounds(:)
863
864CALL stat_bin(sample, bin, nbin, start, finish, mask, binbounds)
865mode = rmiss
866IF (ALLOCATED(bin)) THEN
867 loc = maxloc(bin)
868 mode = (binbounds(loc(1)) + binbounds(loc(1)+1))*0.5
869ENDIF
870
871END FUNCTION stat_mode_histogramr
872
873
874FUNCTION stat_mode_histogramd(sample, nbin, start, finish, mask) RESULT(mode)
875DOUBLE PRECISION,INTENT(in) :: sample(:)
876INTEGER,INTENT(in) :: nbin
877DOUBLE PRECISION,INTENT(in),OPTIONAL :: start
878DOUBLE PRECISION,INTENT(in),OPTIONAL :: finish
879LOGICAL,INTENT(in),OPTIONAL :: mask(:)
880
881DOUBLE PRECISION :: mode
882
883INTEGER :: loc(1)
884INTEGER,ALLOCATABLE :: bin(:)
885DOUBLE PRECISION,ALLOCATABLE :: binbounds(:)
886
887CALL stat_bin(sample, bin, nbin, start, finish, mask, binbounds)
888mode = rmiss
889IF (ALLOCATED(bin)) THEN
890 loc = maxloc(bin)
891 mode = (binbounds(loc(1)) + binbounds(loc(1)+1))*0.5d0
892ENDIF
893
894END FUNCTION stat_mode_histogramd
895
896!!$ Calcolo degli NDI calcolati
897!!$
898!!$ Gli NDI vengono calcolati all’interno degli intervalli dei percentili.
899!!$ Questo significa che il numero di NDI è pari al numero di percentili
900!!$ non degeneri meno 1; il metodo di calcolo è il seguente:
901!!$
902!!$ si calcola il Density Index (DI) per ogni intervallo di percentili
903!!$ come rapporto fra quanti intervalli di percentile vengono utilizzati
904!!$ per definire la larghezza dell’intervallo, più di uno se si hanno
905!!$ percentili degeneri, e la larghezza dell’intervallo stesso.
906!!$
907!!$ Si calcola il valore dell’ADI (Actual Density Index) dividendo i DI
908!!$ ottenuti in intervalli aventi percentili degeneri per il numero di
909!!$ intervalli che hanno contribuito a definire quell’intervallo.
910!!$
911!!$ si calcola il valore del 50-esimo percentile dei valori degli ADI
912!!$ precedentemente ottenuti;
913!!$
914!!$ si divide il valore degli ADI per il 50-percentile degli ADI in
915!!$ maniera tale da normalizzarli ottenendo così gli NADI (Normalized
916!!$ Actual Density Index) calcolati
917!!$
918!!$ Si sommano i NADI relativi ad intervalli di percentili degeneri
919!!$ ottenendo gli NDI
920
921
922
923!!$! Sobroutine calculating the number of NDI values
924!!$SUBROUTINE NOFNDI_old(npclint, pcl0, pclint, pcl100, nndi)
925!!$
926!!$INTEGER, INTENT(IN) :: npclint !< number of pclint
927!!$REAL, INTENT(IN) :: pcl0 !< the 0-th percentile = min(sample)
928!!$REAL, DIMENSION(npclint), INTENT(IN) :: pclint !< the percentiles between the 0-th percentile and the 100-th percentile
929!!$REAL, INTENT(IN) :: pcl100 !< the 100-th percentile = max(sample)
930!!$INTEGER, INTENT(OUT) :: nndi !< number of NDI
931!!$
932!!$REAL, DIMENSION(:), allocatable :: pcl !< array inglobing in itself pcl0, pclint and pcl100
933!!$INTEGER :: &
934!!$ npcl,&
935!!$ ipclint,&! integer loop variable on npclint
936!!$ ipcl, & ! integer loop variable on npcl=2+npclint
937!!$ npcl_act ! number of non redundant values in pcl equal to the number of
938!!$ ! pcl_act values
939!!$
940!!$! Calculation of npcl
941!!$npcl=SIZE(pclint)+2
942!!$
943!!$! Allocation of pcl and its initialisation
944!!$ALLOCATE(pcl(npcl))
945!!$pcl(1)=pcl0
946!!$DO ipclint=1,npclint
947!!$ pcl(ipclint+1)=pclint(ipclint)
948!!$ENDDO
949!!$pcl(SIZE(pclint)+2)=pcl100
950!!$
951!!$IF ( ANY(pcl == rmiss) ) THEN
952!!$
953!!$ nndi=0
954!!$
955!!$ELSE
956!!$
957!!$ ! Calculation of npcl_act
958!!$ npcl_act=1
959!!$ DO ipcl=2,npcl
960!!$ IF (pcl(ipcl) /= pcl(ipcl-1)) THEN
961!!$ npcl_act=npcl_act+1
962!!$ ENDIF
963!!$ ENDDO
964!!$
965!!$ ! Definition of number of ndi that is the wanted output value
966!!$ nndi=npcl_act-1
967!!$
968!!$ENDIF
969!!$
970!!$
971!!$END SUBROUTINE NOFNDI_old
972!!$
973!!$
974!!$
975!!$! Sobroutine calculating the ndi values
976!!$SUBROUTINE NDIC_old(npclint, pcl0, pclint, pcl100, nndi, ndi, limbins)
977!!$
978!!$INTEGER, INTENT(IN) :: npclint !< number of pclint
979!!$REAL, INTENT(IN) :: pcl0 !< the 0-th percentile = min(sample)
980!!$REAL, DIMENSION(npclint), INTENT(IN) :: pclint !< the percentiles between the 0-th percentile and the 100-th percentile
981!!$REAL, INTENT(IN) :: pcl100 !< the 100-th percentile = max(sample)
982!!$INTEGER, INTENT(IN) :: nndi !< number of ndi
983!!$REAL, DIMENSION(nndi), INTENT(OUT) :: ndi !< ndi that I want to calculate
984!!$REAL, DIMENSION(nndi+1), INTENT(OUT) :: limbins !< ndi that I want to calculate
985!!$
986!!$REAL, DIMENSION(:), allocatable :: &
987!!$ pcl, & ! array inglobing in itself pcl0, pclint and pcl100
988!!$ pcl_act, & ! actual values of pcls removing redundat information
989!!$ ranges, & ! distances between the different pcl_act(s)
990!!$ di, & ! density indexes
991!!$ adi, & ! actual density indexes
992!!$ nadi ! normalized actual density indexes
993!!$
994!!$INTEGER, DIMENSION(:), allocatable :: weights ! weights = number of redundance times of a certain values
995!!$ ! in pcl values
996!!$
997!!$REAL, DIMENSION(1) :: &
998!!$ perc_vals, & ! perc_vals contains the percentile position (0-100)
999!!$ med ! mediana value
1000!!$
1001!!$REAL, DIMENSION(:), allocatable :: infopclval
1002!!$INTEGER, DIMENSION(:), allocatable :: infopclnum
1003!!$
1004!!$INTEGER :: &
1005!!$ nbins,& ! number of intervals
1006!!$ ibins,& ! integer loop variable on nbins_act
1007!!$ nbins_act,& ! number of non redundant intervals
1008!!$ ibins_act,& ! integer loop variable on nbins_act
1009!!$ ipclint,& ! integer loop variable on npclint
1010!!$ npcl, & ! number of percentiles
1011!!$ ipcl, & ! integer loop variable on npcl
1012!!$ npcl_act,& ! number of non redundant percentiles
1013!!$ ipcl_act,& ! integer loop variable on npcl_act
1014!!$ npclplus, & ! plus number information
1015!!$ ncount ! number of redundant percentiles values
1016!!$
1017!!$
1018!!$! Actual number of percentiles
1019!!$npcl=SIZE(pclint)+2
1020!!$
1021!!$! Allocation of infopclval and infopclnum
1022!!$ALLOCATE(infopclval(npcl))
1023!!$ALLOCATE(infopclnum(npcl))
1024!!$infopclval(:)=0
1025!!$infopclnum(:)=0
1026!!$
1027!!$! Allocation of pcl
1028!!$ALLOCATE(pcl(npcl))
1029!!$! and storing of pcl values
1030!!$pcl(1)=pcl0
1031!!$DO ipclint=1,npclint
1032!!$ pcl(ipclint+1)=pclint(ipclint)
1033!!$ENDDO
1034!!$pcl(SIZE(pclint)+2)=pcl100
1035!!$
1036!!$ ! Calculation of non redundant values of percentiles
1037!!$
1038!!$npcl_act=1
1039!!$infopclval(1) = pcl(1)
1040!!$infopclnum(1) = 1
1041!!$
1042!!$DO ipcl=2,npcl
1043!!$ infopclval(ipcl) = pcl(ipcl)
1044!!$ IF ( pcl(ipcl) /= pcl(ipcl-1) ) THEN
1045!!$ npcl_act = npcl_act + 1
1046!!$ ENDIF
1047!!$ infopclnum(ipcl) = npcl_act
1048!!$ENDDO
1049!!$
1050!!$ ! Allocation of pcl_act
1051!!$ALLOCATE(pcl_act(npcl_act))
1052!!$ ! and storing in pcl_act of percentiles values
1053!!$DO ipcl_act=1,npcl_act
1054!!$ DO ipcl=1,npcl
1055!!$ IF (ipcl_act == infopclnum(ipcl)) THEN
1056!!$ pcl_act(ipcl_act) = pcl(ipcl)
1057!!$ CYCLE
1058!!$ ENDIF
1059!!$ ENDDO
1060!!$ENDDO
1061!!$
1062!!$
1063!!$ ! Allocation of ranges, weights and di and their initialisation
1064!!$ALLOCATE(ranges(npcl_act-1))
1065!!$ALLOCATE(weights(npcl_act-1))
1066!!$ALLOCATE(di(npcl_act-1))
1067!!$ranges(:)=0
1068!!$di(:)=0
1069!!$weights(:)=0
1070!!$
1071!!$ ! Definition of nbins_act
1072!!$nbins_act=npcl_act-1
1073!!$ ! Cycle on ibins_act for calculating ranges and weights
1074!!$ ! and consequently the di values
1075!!$DO ibins_act=1,nbins_act
1076!!$ ranges(ibins_act)=pcl_act(ibins_act+1) - pcl_act(ibins_act)
1077!!$
1078!!$ IF ( pcl_act(ibins_act+1) == pcl_act(npcl_act) ) THEN
1079!!$ weights(ibins_act) = COUNT( ibins_act == infopclnum ) + &
1080!!$ COUNT( ibins_act+1 == infopclnum ) - 1
1081!!$ ELSE
1082!!$ weights(ibins_act) = COUNT ( ibins_act == infopclnum )
1083!!$ ENDIF
1084!!$ di(ibins_act) = weights(ibins_act)/ranges(ibins_act)
1085!!$ENDDO
1086!!$
1087!!$ ! Allocation of adi and its initialisation
1088!!$ALLOCATE(adi(npcl-1))
1089!!$adi(:)=0
1090!!$npclplus=0
1091!!$DO ibins_act=1,nbins_act
1092!!$ ncount=weights(ibins_act)
1093!!$ ! Calculation of adi
1094!!$ DO ibins=npclplus + 1,npclplus + ncount
1095!!$ adi(ibins)=di(ibins_act)/ncount
1096!!$ ENDDO
1097!!$ npclplus=npclplus+ncount
1098!!$ENDDO
1099!!$
1100!!$ ! Mediana calculation for perc_vals
1101!!$perc_vals(1)=50
1102!!$med = stat_percentile(sample=adi, perc_vals=perc_vals)
1103!!$ ! Allocation of nadi and its initialisation
1104!!$ALLOCATE(nadi(npcl-1))
1105!!$nadi(:)=0
1106!!$ ! Definition of values of nadi
1107!!$nadi(:) = adi(:)/med(1)
1108!!$
1109!!$ ! Initialisation of ndi
1110!!$ndi(:)=0
1111!!$ ! Calculation of the ndi values
1112!!$ipcl_act=1
1113!!$nbins=npcl-1
1114!!$DO ibins=1,nbins
1115!!$ ndi(ipcl_act)=nadi(ibins)+ndi(ipcl_act)
1116!!$ IF ( ( pcl(ibins+1) /= pcl(ibins) ) .AND. &
1117!!$ ( pcl(ibins+1) /= pcl(npcl) ) ) THEN
1118!!$ ipcl_act=ipcl_act+1
1119!!$ ENDIF
1120!!$ENDDO
1121!!$
1122!!$DO ipcl_act=1,npcl_act
1123!!$ limbins(ipcl_act)=pcl_act(ipcl_act)
1124!!$ENDDO
1125!!$
1126!!$ ! Deallocation part
1127!!$DEALLOCATE(infopclval)
1128!!$DEALLOCATE(infopclnum)
1129!!$DEALLOCATE(pcl)
1130!!$DEALLOCATE(pcl_act)
1131!!$DEALLOCATE(ranges)
1132!!$DEALLOCATE(weights)
1133!!$DEALLOCATE(di)
1134!!$DEALLOCATE(adi)
1135!!$DEALLOCATE(nadi)
1136!!$
1137!!$END SUBROUTINE NDIC_old
1138!!$
1139!!$
1140!!$!> Calculate the number of NDI values
1141!!$SUBROUTINE NOFNDI(pcl, nndi)
1142!!$
1143!!$REAL, DIMENSION(:), INTENT(IN) :: pcl !< the percentiles between the 0-th percentile and the 100-th percentile
1144!!$INTEGER, INTENT(OUT) :: nndi !< number of NDI
1145!!$
1146!!$IF ( ANY(pcl == rmiss) ) then
1147!!$ nndi=0
1148!!$else
1149!!$ nndi=count_distinct(pcl)-1
1150!!$end IF
1151!!$
1152!!$END SUBROUTINE NOFNDI
1153
1154
1155
1156!!$!example to manage exceptions
1157!!$
1158!!$use,intrinsic :: IEEE_EXCEPTIONS
1159!!$
1160!!$logical fail_o,fail_zero
1161!!$
1162!!$a=10.
1163!!$b=tiny(0.)
1164!!$
1165!!$call safe_divide(a,b,c,fail_o,fail_zero)
1166!!$
1167!!$print*,fail_o,fail_zero
1168!!$print *,a,b,c
1169!!$
1170!!$b=0.
1171!!$call safe_divide(a,b,c,fail_o,fail_zero)
1172!!$
1173!!$print*,fail_o,fail_zero
1174!!$print *,a,b,c
1175!!$
1176!!$contains
1177!!$
1178!!$subroutine safe_divide(a, b, c, fail_o,fail_zero)
1179!!$
1180!!$real a, b, c
1181!!$logical fail_o,fail_zero
1182!!$type(IEEE_STATUS_TYPE) status
1183!!$! save the current floating-point environment, turn halting for
1184!!$! divide-by-zero off, and clear any previous divide-by-zero flag
1185!!$call IEEE_GET_STATUS(status)
1186!!$call IEEE_SET_HALTING_MODE(IEEE_DIVIDE_BY_ZERO, .false.)
1187!!$call IEEE_SET_HALTING_MODE(ieee_overflow, .false.)
1188!!$
1189!!$call IEEE_SET_FLAG(ieee_overflow, .false.)
1190!!$call IEEE_SET_FLAG(IEEE_DIVIDE_BY_ZERO, .false.)
1191!!$! perform the operation
1192!!$c = a/b
1193!!$! determine if a failure occurred and restore the floating-point environment
1194!!$call IEEE_GET_FLAG(ieee_overflow, fail_o)
1195!!$call IEEE_GET_FLAG(IEEE_DIVIDE_BY_ZERO, fail_zero)
1196!!$call IEEE_SET_STATUS(status)
1197!!$end subroutine safe_divide
1198!!$end program
1199!!$
1200
1201
1202! here you have to adopt the example above in the code below the use the wrong logic isnan()
1203
1204!!$!check NaN
1205!!$if (any(isnan(di)) .and. .not. all(isnan(di))) then
1206!!$
1207!!$! from left
1208!!$ do i=2,size(delta)
1209!!$ if (isnan(di(i)) .and. .not. isnan(di(i-1))) then
1210!!$ delta(i) = delta(i-1)
1211!!$ w(i) = w(i) + w(i-1)
1212!!$ end if
1213!!$ end do
1214!!$
1215!!$! recompute
1216!!$ print *,"WW=",w
1217!!$ di = w/delta
1218!!$
1219!!$ if (any(isnan(di)) .and. .not. all(isnan(di))) then
1220!!$
1221!!$! from right
1222!!$ do i=size(delta)-1,1,-1
1223!!$ if (isnan(di(i)) .and. .not. isnan(di(i+1))) then
1224!!$ delta(i) = delta(i+1)
1225!!$ w(i) = w(i) + w(i+1)
1226!!$ end if
1227!!$ end do
1228!!$
1229!!$! one more step
1230!!$ call DensityIndex(di,perc_vals,limbins)
1231!!$ end if
1232!!$end if
1233
1234
1235!!$subroutine DensityIndex_old(di,perc_vals,limbins)
1236!!$real,intent(inout) :: di(:)
1237!!$real,intent(in) :: perc_vals(:)
1238!!$real,intent(in) :: limbins(:)
1239!!$
1240!!$real :: delta(size(di)),w(size(di))
1241!!$integer :: i
1242!!$
1243!!$do i=1,size(delta)
1244!!$ delta(i) = limbins(i+1) - limbins(i)
1245!!$ w(i) = perc_vals(i+1) - perc_vals(i)
1246!!$end do
1247!!$
1248!!$di=rmiss
1249!!$
1250!!$if ( .not. all(delta == 0.)) then
1251!!$ call DensityIndex_recurse(delta,w)
1252!!$ di = w/delta
1253!!$end if
1254!!$
1255!!$end subroutine DensityIndex_old
1256!!$
1257!!$recursive subroutine DensityIndex_recurse(delta,w)
1258!!$real :: delta(:),w(:)
1259!!$integer :: i
1260!!$
1261!!$!check divide by 0
1262!!$if (any(delta == 0.0)) then
1263!!$
1264!!$! from left
1265!!$ do i=2,size(delta)
1266!!$ if (delta(i) == 0.0 .and. delta(i-1) > 0.0 ) then
1267!!$ delta(i) = delta(i-1)
1268!!$ w(i) = w(i) + w(i-1)
1269!!$ end if
1270!!$ end do
1271!!$
1272!!$! from right
1273!!$ do i=size(delta)-1,1,-1
1274!!$ if (delta(i) == 0.0 .and. delta(i+1) > 0.0 )then
1275!!$ delta(i) = delta(i+1)
1276!!$ w(i) = w(i) + w(i+1)
1277!!$ end if
1278!!$ end do
1279!!$
1280!!$end if
1281!!$
1282!!$! one more step
1283!!$if (any(delta == 0.0)) then
1284!!$ call DensityIndex_recurse(delta,w)
1285!!$end if
1286!!$
1287!!$end subroutine DensityIndex_recurse
1288!!$
1289
1290subroutine densityindex(di,nlimbins,occu,rnum,limbins)
1291real,intent(out) :: di(:)
1292real,intent(out) :: nlimbins(:)
1293integer,intent(out) :: occu(:)
1294REAL, DIMENSION(:), INTENT(IN) :: rnum !< data to analize
1295real,intent(in) :: limbins(:)
1296
1297real :: nnum(size(rnum))
1298integer :: i,k,sample_count
1299logical :: sample_mask(size(rnum))
1300
1301di=rmiss
1302occu=imiss
1303nlimbins=rmiss
1304
1305nlimbins(1)=limbins(1) ! compute unique limits
1306k=1
1307do i=2,size(limbins)
1308 if (limbins(i) /= limbins(k)) then
1309 k=k+1
1310 nlimbins(k)= limbins(i)
1311 end if
1312end do
1313
1314di=rmiss
1315if (k == 1) return
1316
1317sample_mask = (rnum /= rmiss) ! remove missing values
1318sample_count = count(sample_mask)
1319IF (sample_count == 0) RETURN
1320nnum(1:sample_count) = pack(rnum, mask=sample_mask)
1321
1322do i=1,k-2 ! compute occorrence and density index
1323 occu(i)=count(nnum>=nlimbins(i) .and. nnum<nlimbins(i+1))
1324 di(i) = float(occu(i)) / (nlimbins(i+1) - nlimbins(i))
1325end do
1326
1327i=k-1 ! the last if is <=
1328occu(i)=count(nnum>=nlimbins(i) .and. nnum<=nlimbins(i+1))
1329di(i) = float(occu(i)) / (nlimbins(i+1) - nlimbins(i))
1330
1331end subroutine densityindex
1332
1333
1334!> Compute Normalized Density Index
1335SUBROUTINE normalizeddensityindex(rnum, perc_vals, ndi, nlimbins)
1336
1337REAL, DIMENSION(:), INTENT(IN) :: rnum !< data to analize
1338REAL, DIMENSION(:), INTENT(IN) :: perc_vals !<the percentiles values to be computed, between 0. and 100.
1339REAL, DIMENSION(:), INTENT(OUT) :: ndi !< normalized density index
1340REAL, DIMENSION(:), INTENT(OUT) :: nlimbins !< the extreme values of data taken in account for ndi computation
1341
1342REAL, DIMENSION(size(ndi)) :: di
1343INTEGER, DIMENSION(size(ndi)) :: occu
1344REAL, DIMENSION(size(nlimbins)) :: limbins
1345real :: med
1346integer :: i,k,middle
1347
1348ndi=rmiss
1349limbins = stat_percentile(rnum,perc_vals) ! compute percentile
1350
1351call densityindex(di,nlimbins,occu,rnum,limbins)
1352
1353! Mediana calculation for density index
1354k=0
1355middle=count(c_e(rnum))/2
1356
1357do i=1,count(c_e(occu))
1358 k=k+occu(i)
1359 if (k > middle) then
1360 if (k > 1 .and. (k - occu(i)) == middle) then
1361 med = (di(i-1) + di(i)) / 2.
1362 else
1363 med = di(i)
1364 end if
1365 exit
1366 end if
1367end do
1368
1369!weighted density index
1370ndi(:count(c_e(di))) = min(pack(di,mask=c_e(di))/med,1.0)
1371
1372END SUBROUTINE normalizeddensityindex
1373
1374
1375!! TODO translate from python
1376
1377!!$def gauss(x, A=1, mu=1, sigma=1):
1378!!$ """
1379!!$ Evaluate Gaussian.
1380!!$
1381!!$ Parameters
1382!!$ ----------
1383!!$ A : float
1384!!$ Amplitude.
1385!!$ mu : float
1386!!$ Mean.
1387!!$ std : float
1388!!$ Standard deviation.
1389!!$
1390!!$ """
1391!!$ return np.real(A * np.exp(-(x - mu)**2 / (2 * sigma**2)))
1392!!$
1393!!$def fit_direct(x, y, F=0, weighted=True, _weights=None):
1394!!$ """Fit a Gaussian to the given data.
1395!!$
1396!!$ Returns a fit so that y ~ gauss(x, A, mu, sigma)
1397!!$
1398!!$ Parameters
1399!!$ ----------
1400!!$ x : ndarray
1401!!$ Sampling positions.
1402!!$ y : ndarray
1403!!$ Sampled values.
1404!!$ F : float
1405!!$ Ignore values of y <= F.
1406!!$ weighted : bool
1407!!$ Whether to use weighted least squares. If True, weigh
1408!!$ the error function by y, ensuring that small values
1409!!$ has less influence on the outcome.
1410!!$
1411!!$ Additional Parameters
1412!!$ ---------------------
1413!!$ _weights : ndarray
1414!!$ Weights used in weighted least squares. For internal use
1415!!$ by fit_iterative.
1416!!$
1417!!$ Returns
1418!!$ -------
1419!!$ A : float
1420!!$ Amplitude.
1421!!$ mu : float
1422!!$ Mean.
1423!!$ std : float
1424!!$ Standard deviation.
1425!!$
1426!!$ """
1427!!$ mask = (y > F)
1428!!$ x = x[mask]
1429!!$ y = y[mask]
1430!!$
1431!!$ if _weights is None:
1432!!$ _weights = y
1433!!$ else:
1434!!$ _weights = _weights[mask]
1435!!$
1436!!$ # We do not want to risk working with negative values
1437!!$ np.clip(y, 1e-10, np.inf, y)
1438!!$
1439!!$ e = np.ones(len(x))
1440!!$ if weighted:
1441!!$ e = e * (_weights**2)
1442!!$
1443!!$ v = (np.sum(np.vander(x, 5) * e[:, None], axis=0))[::-1]
1444!!$ A = v[sl.hankel([0, 1, 2], [2, 3, 4])]
1445!!$
1446!!$ ly = e * np.log(y)
1447!!$ ls = np.sum(ly)
1448!!$ x_ls = np.sum(ly * x)
1449!!$ xx_ls = np.sum(ly * x**2)
1450!!$ B = np.array([ls, x_ls, xx_ls])
1451!!$
1452!!$ (a, b, c), res, rank, s = np.linalg.lstsq(A, B)
1453!!$
1454!!$ A = np.exp(a - (b**2 / (4 * c)))
1455!!$ mu = -b / (2 * c)
1456!!$ sigma = sp.sqrt(-1 / (2 * c))
1457!!$
1458!!$ return A, mu, sigma
1459!!$
1460!!$def fit_iterative(x, y, F=0, weighted=True, N=10):
1461!!$ """Fit a Gaussian to the given data.
1462!!$
1463!!$ Returns a fit so that y ~ gauss(x, A, mu, sigma)
1464!!$
1465!!$ This function iteratively fits using fit_direct.
1466!!$
1467!!$ Parameters
1468!!$ ----------
1469!!$ x : ndarray
1470!!$ Sampling positions.
1471!!$ y : ndarray
1472!!$ Sampled values.
1473!!$ F : float
1474!!$ Ignore values of y <= F.
1475!!$ weighted : bool
1476!!$ Whether to use weighted least squares. If True, weigh
1477!!$ the error function by y, ensuring that small values
1478!!$ has less influence on the outcome.
1479!!$ N : int
1480!!$ Number of iterations.
1481!!$
1482!!$ Returns
1483!!$ -------
1484!!$ A : float
1485!!$ Amplitude.
1486!!$ mu : float
1487!!$ Mean.
1488!!$ std : float
1489!!$ Standard deviation.
1490!!$
1491!!$ """
1492!!$ y_ = y
1493!!$ for i in range(N):
1494!!$ p = fit_direct(x, y, weighted=True, _weights=y_)
1495!!$ A, mu, sigma = p
1496!!$ y_ = gauss(x, A, mu, sigma)
1497!!$
1498!!$ return np.real(A), np.real(mu), np.real(sigma)
1499!!$
1500
1501
1502
1503END MODULE simple_stat
Function to check whether a value is missing or not.
Compute the average of the random variable provided, taking into account missing data.
Bin a sample into equally spaced intervals to form a histogram.
Compute the linear correlation coefficient between the two random variables provided,...
Compute the linear regression coefficients between the two random variables provided,...
Compute the mode of the random variable provided taking into account missing data.
Compute a set of percentiles for a random variable.
Compute the standard deviation of the random variable provided, taking into account missing data.
Compute the variance of the random variable provided, taking into account missing data.
This module defines usefull general purpose function and subroutine.
Definitions of constants and functions for working with missing values.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Module for basic statistical computations taking into account missing data.

Generated with Doxygen.