40 MODULE PROCEDURE stat_averager, stat_averaged
55 MODULE PROCEDURE stat_variancer, stat_varianced
70 MODULE PROCEDURE stat_stddevr, stat_stddevd
90 MODULE PROCEDURE stat_linear_corrr, stat_linear_corrd
108 MODULE PROCEDURE stat_linear_regressionr, stat_linear_regressiond
122 MODULE PROCEDURE stat_percentiler, stat_percentiled
142 MODULE PROCEDURE stat_binr, stat_bind
161 MODULE PROCEDURE stat_mode_histogramr, stat_mode_histogramd
167 normalizeddensityindex
172FUNCTION stat_averager(sample, mask, nomiss)
RESULT(average)
173REAL,
INTENT(in) :: sample(:)
174LOGICAL,
OPTIONAL,
INTENT(in) :: mask(:)
175LOGICAL,
OPTIONAL,
INTENT(in) :: nomiss
179INTEGER :: sample_count
180LOGICAL :: sample_mask(SIZE(sample))
182IF (optio_log(nomiss))
THEN
183 average = sum(sample)/
SIZE(sample)
185 sample_mask = (sample /= rmiss)
186 IF (
PRESENT(mask)) sample_mask = sample_mask .AND. mask
187 sample_count = count(sample_mask)
189 IF (sample_count > 0)
THEN
191 average = sum(sample, mask=sample_mask)/sample_count
197END FUNCTION stat_averager
200FUNCTION stat_averaged(sample, mask, nomiss)
RESULT(average)
201DOUBLE PRECISION,
INTENT(in) :: sample(:)
202LOGICAL,
OPTIONAL,
INTENT(in) :: mask(:)
203LOGICAL,
OPTIONAL,
INTENT(in) :: nomiss
205DOUBLE PRECISION :: average
207INTEGER :: sample_count
208LOGICAL :: sample_mask(SIZE(sample))
210IF (optio_log(nomiss))
THEN
211 average = sum(sample)/
SIZE(sample)
213 sample_mask = (sample /= dmiss)
214 IF (
PRESENT(mask)) sample_mask = sample_mask .AND. mask
215 sample_count = count(sample_mask)
217 IF (sample_count > 0)
THEN
219 average = sum(sample, mask=sample_mask)/sample_count
225END FUNCTION stat_averaged
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
238INTEGER :: sample_count, i
239LOGICAL :: sample_mask(SIZE(sample))
241IF (optio_log(nomiss))
THEN
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)
248 variance = sum((sample-laverage)**2)/
SIZE(sample)
252 sample_mask = (sample /= rmiss)
253 IF (
PRESENT(mask)) sample_mask = sample_mask .AND. mask
254 sample_count = count(sample_mask)
256 IF (sample_count > 0)
THEN
258 laverage = sum(sample, mask=sample_mask)/sample_count
259 IF (
PRESENT(average)) average = laverage
262 DO i = 1,
SIZE(sample)
263 IF (sample_mask(i)) variance = variance + (sample(i)-laverage)**2
265 IF (optio_log(nm1))
THEN
266 variance = variance/max(sample_count-1,1)
268 variance = variance/sample_count
271 IF (
PRESENT(average)) average = rmiss
277END FUNCTION stat_variancer
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
287DOUBLE PRECISION :: variance
289DOUBLE PRECISION :: laverage
290INTEGER :: sample_count, i
291LOGICAL :: sample_mask(SIZE(sample))
293IF (optio_log(nomiss))
THEN
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)
300 variance = sum((sample-laverage)**2)/
SIZE(sample)
304 sample_mask = (sample /= dmiss)
305 IF (
PRESENT(mask)) sample_mask = sample_mask .AND. mask
306 sample_count = count(sample_mask)
308 IF (sample_count > 0)
THEN
310 laverage = sum(sample, mask=sample_mask)/sample_count
311 IF (
PRESENT(average)) average = laverage
314 DO i = 1,
SIZE(sample)
315 IF (sample_mask(i)) variance = variance + (sample(i)-laverage)**2
317 IF (optio_log(nm1))
THEN
318 variance = variance/max(sample_count-1,1)
320 variance = variance/sample_count
323 IF (
PRESENT(average)) average = dmiss
329END FUNCTION stat_varianced
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
342IF (
c_e(stddev)) stddev = sqrt(stddev)
344END FUNCTION stat_stddevr
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
354DOUBLE PRECISION :: stddev
357IF (
c_e(stddev)) stddev = sqrt(stddev)
359END FUNCTION stat_stddevd
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
375REAL :: laverage1, laverage2, lvariance1, lvariance2
376INTEGER :: sample_count, i
377LOGICAL :: sample_mask(SIZE(sample1))
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
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
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
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
406 lvariance1 = lvariance1/sample_count
407 lvariance2 = lvariance2/sample_count
408 IF (
PRESENT(variance1)) variance1 = lvariance1
409 IF (
PRESENT(variance2)) variance2 = lvariance2
412 DO i = 1,
SIZE(sample1)
413 IF (sample_mask(i)) linear_corr = linear_corr + &
414 (sample1(i)-laverage1)*(sample2(i)-laverage2)
416 linear_corr = linear_corr/sample_count / sqrt(lvariance1*lvariance2)
418 IF (
PRESENT(average1)) average1 = rmiss
419 IF (
PRESENT(average2)) average2 = rmiss
420 IF (
PRESENT(variance1)) variance1 = rmiss
421 IF (
PRESENT(variance2)) variance2 = rmiss
425END FUNCTION stat_linear_corrr
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
439DOUBLE PRECISION :: linear_corr
441DOUBLE PRECISION :: laverage1, laverage2, lvariance1, lvariance2
442INTEGER :: sample_count, i
443LOGICAL :: sample_mask(SIZE(sample1))
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
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
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
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
472 lvariance1 = lvariance1/sample_count
473 lvariance2 = lvariance2/sample_count
474 IF (
PRESENT(variance1)) variance1 = lvariance1
475 IF (
PRESENT(variance2)) variance2 = lvariance2
478 DO i = 1,
SIZE(sample1)
479 IF (sample_mask(i)) linear_corr = linear_corr + &
480 (sample1(i)-laverage1)*(sample2(i)-laverage2)
482 linear_corr = linear_corr/sample_count / sqrt(lvariance1*lvariance2)
484 IF (
PRESENT(average1)) average1 = dmiss
485 IF (
PRESENT(average2)) average2 = dmiss
486 IF (
PRESENT(variance1)) variance1 = dmiss
487 IF (
PRESENT(variance2)) variance2 = dmiss
491END FUNCTION stat_linear_corrd
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(:)
501REAL :: laverage1, laverage2
502INTEGER :: sample_count
503LOGICAL :: sample_mask(SIZE(sample1))
505IF (
SIZE(sample1) /=
SIZE(sample2))
THEN
511sample_mask = (sample1 /= rmiss .AND. sample2 /= rmiss)
512IF (
PRESENT(mask)) sample_mask = sample_mask .AND. mask
513sample_count = count(sample_mask)
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
527END SUBROUTINE stat_linear_regressionr
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(:)
537DOUBLE PRECISION :: laverage1, laverage2
538INTEGER :: sample_count
539LOGICAL :: sample_mask(SIZE(sample1))
541IF (
SIZE(sample1) /=
SIZE(sample2))
THEN
547sample_mask = (sample1 /= dmiss .AND. sample2 /= dmiss)
548IF (
PRESENT(mask)) sample_mask = sample_mask .AND. mask
549sample_count = count(sample_mask)
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
563END SUBROUTINE stat_linear_regressiond
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))
573REAL :: lsample(SIZE(sample)), rindex
574INTEGER :: sample_count, j, iindex
575LOGICAL :: sample_mask(SIZE(sample))
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
584 sample_count =
SIZE(sample)
587IF (sample_count ==
SIZE(sample))
THEN
590 lsample(1:sample_count) = pack(sample, mask=sample_mask)
593IF (sample_count == 1)
THEN
594 percentile(:) = lsample(1)
600CALL sort(lsample(1:sample_count))
612DO j = 1,
SIZE(perc_vals)
613 IF (perc_vals(j) >= 0. .AND. perc_vals(j) <= 100.)
THEN
615 rindex = real(sample_count-1, kind=kind(rindex))*perc_vals(j)/100.+1.
617 iindex = min(max(int(rindex), 1), sample_count-1)
620 percentile(j) = lsample(iindex)*(real(iindex+1, kind=kind(rindex))-rindex) &
621 + lsample(iindex+1)*(rindex-real(iindex, kind=kind(rindex)))
626END FUNCTION stat_percentiler
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))
636DOUBLE PRECISION :: lsample(SIZE(sample)), rindex
637INTEGER :: sample_count, j, iindex
638LOGICAL :: sample_mask(SIZE(sample))
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
648 sample_count =
SIZE(sample)
651IF (sample_count ==
SIZE(sample))
THEN
654 lsample(1:sample_count) = pack(sample, mask=sample_mask)
657IF (sample_count == 1)
THEN
658 percentile(:) = lsample(1)
663CALL sort(lsample(1:sample_count))
675DO j = 1,
SIZE(perc_vals)
676 IF (perc_vals(j) >= 0.d0 .AND. perc_vals(j) <= 100.d0)
THEN
678 rindex = real(sample_count-1, kind=kind(rindex))*perc_vals(j)/100.d0+1.d0
680 iindex = min(max(int(rindex), 1), sample_count-1)
682 percentile(j) = lsample(iindex)*(real(iindex+1, kind=kind(rindex))-rindex) &
683 + lsample(iindex+1)*(rindex-real(iindex, kind=kind(rindex)))
687END FUNCTION stat_percentiled
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(:)
700REAL :: lstart, lfinish, incr
701REAL,
ALLOCATABLE :: lbinbounds(:)
702LOGICAL :: lmask(SIZE(sample))
706IF (
PRESENT(mask))
THEN
711lmask = lmask .AND.
c_e(sample)
712IF (count(lmask) < 1)
RETURN
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
720incr = (lfinish-lstart)/nbin
723ALLOCATE(lbinbounds(nbin+1))
726 lbinbounds(i) = lstart + (i-1)*incr
728lbinbounds(nbin+1) = lfinish
731 bin(i) = count(sample >= lbinbounds(i) .AND. sample < lbinbounds(i+1) .AND. lmask)
734bin(nbin) = count(sample >= lbinbounds(nbin) .AND. sample <= lbinbounds(nbin+1) .AND. lmask)
737IF (
PRESENT(binbounds)) binbounds = lbinbounds
739END SUBROUTINE stat_binr
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(:)
753REAL :: lstart, lfinish, incr
754LOGICAL :: lmask(SIZE(sample))
758IF (
PRESENT(mask))
THEN
763lmask = lmask .AND.
c_e(sample)
764IF (count(lmask) < 1)
RETURN
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
772incr = (lfinish-lstart)/nbin
777DO i = 1,
SIZE(sample)
779 ind = int((sample(i)-lstart)/incr) + 1
780 IF (ind > 0 .AND. ind <= nbin)
THEN
781 bin(ind) = bin(ind) + 1
783 IF (sample(i) == finish) bin(nbin) = bin(nbin) + 1
788IF (
PRESENT(binbounds))
THEN
789 ALLOCATE(binbounds(nbin+1))
791 binbounds(i) = lstart + (i-1)*incr
793 binbounds(nbin+1) = lfinish
796END SUBROUTINE stat_binr2
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(:)
809DOUBLE PRECISION :: lstart, lfinish, incr
810DOUBLE PRECISION,
ALLOCATABLE :: lbinbounds(:)
811LOGICAL :: lmask(SIZE(sample))
815IF (
PRESENT(mask))
THEN
820lmask = lmask .AND.
c_e(sample)
821IF (count(lmask) < 1)
RETURN
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
829incr = (lfinish-lstart)/nbin
832ALLOCATE(lbinbounds(nbin+1))
835 lbinbounds(i) = lstart + (i-1)*incr
837lbinbounds(nbin+1) = lfinish
840 bin(i) = count(sample >= lbinbounds(i) .AND. sample < lbinbounds(i+1) .AND. lmask)
843bin(nbin) = count(sample >= lbinbounds(nbin) .AND. sample <= lbinbounds(nbin+1) .AND. lmask)
846IF (
PRESENT(binbounds)) binbounds = lbinbounds
848END SUBROUTINE stat_bind
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(:)
861INTEGER,
ALLOCATABLE :: bin(:)
862REAL,
ALLOCATABLE :: binbounds(:)
864CALL stat_bin(sample, bin, nbin, start, finish, mask, binbounds)
866IF (
ALLOCATED(bin))
THEN
868 mode = (binbounds(loc(1)) + binbounds(loc(1)+1))*0.5
871END FUNCTION stat_mode_histogramr
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(:)
881DOUBLE PRECISION :: mode
884INTEGER,
ALLOCATABLE :: bin(:)
885DOUBLE PRECISION,
ALLOCATABLE :: binbounds(:)
887CALL stat_bin(sample, bin, nbin, start, finish, mask, binbounds)
889IF (
ALLOCATED(bin))
THEN
891 mode = (binbounds(loc(1)) + binbounds(loc(1)+1))*0.5d0
894END FUNCTION stat_mode_histogramd
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(:)
1297real :: nnum(size(rnum))
1298integer :: i,k,sample_count
1299logical :: sample_mask(size(rnum))
1305nlimbins(1)=limbins(1)
1308 if (limbins(i) /= limbins(k))
then
1310 nlimbins(k)= limbins(i)
1317sample_mask = (rnum /= rmiss)
1318sample_count = count(sample_mask)
1319IF (sample_count == 0)
RETURN
1320nnum(1:sample_count) = pack(rnum, mask=sample_mask)
1323 occu(i)=count(nnum>=nlimbins(i) .and. nnum<nlimbins(i+1))
1324 di(i) = float(occu(i)) / (nlimbins(i+1) - nlimbins(i))
1328occu(i)=count(nnum>=nlimbins(i) .and. nnum<=nlimbins(i+1))
1329di(i) = float(occu(i)) / (nlimbins(i+1) - nlimbins(i))
1331end subroutine densityindex
1335SUBROUTINE normalizeddensityindex(rnum, perc_vals, ndi, nlimbins)
1337REAL,
DIMENSION(:),
INTENT(IN) :: rnum
1338REAL,
DIMENSION(:),
INTENT(IN) :: perc_vals
1339REAL,
DIMENSION(:),
INTENT(OUT) :: ndi
1340REAL,
DIMENSION(:),
INTENT(OUT) :: nlimbins
1342REAL,
DIMENSION(size(ndi)) :: di
1343INTEGER,
DIMENSION(size(ndi)) :: occu
1344REAL,
DIMENSION(size(nlimbins)) :: limbins
1346integer :: i,k,middle
1351call densityindex(di,nlimbins,occu,rnum,limbins)
1355middle=count(
c_e(rnum))/2
1357do i=1,count(
c_e(occu))
1359 if (k > middle)
then
1360 if (k > 1 .and. (k - occu(i)) == middle)
then
1361 med = (di(i-1) + di(i)) / 2.
1370ndi(:count(
c_e(di))) = min(pack(di,mask=
c_e(di))/med,1.0)
1372END SUBROUTINE normalizeddensityindex
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.