libsim Versione 7.2.0
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/>.
29IMPLICIT NONE
30
39INTERFACE stat_average
40 MODULE PROCEDURE stat_averager, stat_averaged
41END INTERFACE
42
55 MODULE PROCEDURE stat_variancer, stat_varianced
56END INTERFACE
57
69INTERFACE stat_stddev
70 MODULE PROCEDURE stat_stddevr, stat_stddevd
71END INTERFACE
72
90 MODULE PROCEDURE stat_linear_corrr, stat_linear_corrd
91END INTERFACE
92
108 MODULE PROCEDURE stat_linear_regressionr, stat_linear_regressiond
109END INTERFACE
110
122 MODULE PROCEDURE stat_percentiler, stat_percentiled
123END INTERFACE
124
141INTERFACE stat_bin
142 MODULE PROCEDURE stat_binr, stat_bind
143END INTERFACE stat_bin
144
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
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
1335SUBROUTINE normalizeddensityindex(rnum, perc_vals, ndi, nlimbins)
1336
1337REAL, DIMENSION(:), INTENT(IN) :: rnum
1338REAL, DIMENSION(:), INTENT(IN) :: perc_vals
1339REAL, DIMENSION(:), INTENT(OUT) :: ndi
1340REAL, DIMENSION(:), INTENT(OUT) :: nlimbins
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.
Definition: simple_stat.f90:39
Bin a sample into equally spaced intervals to form a histogram.
Compute the linear correlation coefficient between the two random variables provided,...
Definition: simple_stat.f90:89
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.
Definition: simple_stat.f90:69
Compute the variance of the random variable provided, taking into account missing data.
Definition: simple_stat.f90:54
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.
Definition: simple_stat.f90:25

Generated with Doxygen.