The ParaMonte Documentation Website
Current view: top level - kernel - CorrCoef_mod.f90 (source / functions) Hit Total Coverage
Test: ParaMonte 1.5.1 :: MPI Parallel Kernel - Code Coverage Report Lines: 37 44 84.1 %
Date: 2021-01-08 13:07:16 Functions: 2 2 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       2             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       3             : !!!!
       4             : !!!!   MIT License
       5             : !!!!
       6             : !!!!   ParaMonte: plain powerful parallel Monte Carlo library.
       7             : !!!!
       8             : !!!!   Copyright (C) 2012-present, The Computational Data Science Lab
       9             : !!!!
      10             : !!!!   This file is part of the ParaMonte library.
      11             : !!!!
      12             : !!!!   Permission is hereby granted, free of charge, to any person obtaining a 
      13             : !!!!   copy of this software and associated documentation files (the "Software"), 
      14             : !!!!   to deal in the Software without restriction, including without limitation 
      15             : !!!!   the rights to use, copy, modify, merge, publish, distribute, sublicense, 
      16             : !!!!   and/or sell copies of the Software, and to permit persons to whom the 
      17             : !!!!   Software is furnished to do so, subject to the following conditions:
      18             : !!!!
      19             : !!!!   The above copyright notice and this permission notice shall be 
      20             : !!!!   included in all copies or substantial portions of the Software.
      21             : !!!!
      22             : !!!!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 
      23             : !!!!   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 
      24             : !!!!   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. 
      25             : !!!!   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, 
      26             : !!!!   DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR 
      27             : !!!!   OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE 
      28             : !!!!   OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
      29             : !!!!
      30             : !!!!   ACKNOWLEDGMENT
      31             : !!!!
      32             : !!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
      33             : !!!!   As per the ParaMonte library license agreement terms, if you use any parts of 
      34             : !!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your 
      35             : !!!!   work (education/research/industry/development/...) by citing the ParaMonte 
      36             : !!!!   library as described on this page:
      37             : !!!!
      38             : !!!!       https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
      39             : !!!!
      40             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      41             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      42             : 
      43             : !> \brief This module contains procedures for computing various statistical correlation coefficients.
      44             : !> \author Amir Shahmoradi
      45             : 
      46             : module CorrCoef_mod
      47             : 
      48             :     use Constants_mod, only: RK, IK
      49             :     use Err_mod, only: Err_type
      50             :     implicit none
      51             : 
      52             :     character(*), parameter :: MODULE_NAME = "@CorrCoef_mod"
      53             : 
      54             :     type :: CorrCoefSpearman_type
      55             :         integer(IK)             :: ndata
      56             :         real(RK), allocatable   :: Data1(:), Data2(:)
      57             :         real(RK)                :: rho, rhoProb, dStarStar, dStarStarSignif, dStarStarProb
      58             :         type(Err_type)          :: Err
      59             :     contains
      60             :         procedure, nopass :: get => getCorrCoefSpearman
      61             :     end type CorrCoefSpearman_type
      62             : 
      63             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      64             : 
      65             : contains
      66             : 
      67             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      68             : 
      69             :     !> \brief
      70             :     !> Return the Spearman correlation coefficient of the two input data arrays.
      71             :     !> 
      72             :     !> @param[in] ndata             :   The length of the two input data arrays.
      73             :     !> @param[in] Data1             :   The input data array of size `(1:ndata)`.
      74             :     !> @param[in] Data2             :   The input data array of size `(1:ndata)`.
      75             :     !> @param[out] rho              :   The Spearman rank correlation coefficient.
      76             :     !> @param[out] rhoProb          :   The two-sided significance level of the deviation of `rho` from zero.
      77             :     !> @param[out] dStarStar        :   The sum-squared difference of ranks of the two data vectors.
      78             :     !> @param[out] dStarStarSignif  :   The number of standard deviations by which `dStarStar` deviates from its null-hypothesis expected value.
      79             :     !> @param[out] dStarStarProb    :   The two-sided significance level of the deviation represented by `dStarStarSignif`.
      80             :     !> @param[out] Err              :   An object of [Err_type](@ref err_mod::err_type) class, indicating whether an error has occurred.
      81             :     !> 
      82             :     !> \remark
      83             :     !> A small value of either `dStarStarProb` or `rhoProb` indicates a significant correlation.
      84           6 :     subroutine getCorrCoefSpearman(ndata,Data1,Data2,rho,rhoProb,dStarStar,dStarStarSignif,dStarStarProb,Err)
      85             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
      86             :         !DEC$ ATTRIBUTES DLLEXPORT :: getCorrCoefSpearman
      87             : #endif
      88             :         use Constants_mod, only: IK, RK, SQRT2, SPR
      89             :         use Statistics_mod, only: getBetaCDF
      90             :         use Sort_mod, only: sortAscending
      91             : 
      92             :         implicit none
      93             : 
      94             :         integer(IK)     , intent(in)    :: ndata
      95             :         real(RK)        , intent(in)    :: Data1(ndata), Data2(ndata)
      96             :         real(RK)        , intent(out)   :: rho, rhoProb
      97             :         real(RK)        , intent(out)   :: dStarStar, dStarStarSignif, dStarStarProb
      98             :         type(Err_type)  , intent(out)   :: Err
      99             : 
     100             :         character(*)    , parameter     :: PROCEDURE_NAME = MODULE_NAME//"@getCorrCoefSpearman"
     101           6 :         real(RK)                        :: aved,df,en,en3n,fac,sf,sg,t,vard
     102         606 :         real(RK)                        :: WorkSpace1(ndata), WorkSpace2(ndata)
     103             : 
     104         306 :         WorkSpace1 = Data1
     105         306 :         WorkSpace2 = Data2
     106           6 :         call sortAscending(ndata,WorkSpace1,WorkSpace2,Err)
     107           6 :         if (Err%occurred) then
     108             :         ! LCOV_EXCL_START
     109             :             Err%msg = PROCEDURE_NAME//Err%msg
     110             :             return
     111             :         end if
     112             :         ! LCOV_EXCL_STOP
     113           6 :         call crank(ndata,WorkSpace1,sf)
     114             : 
     115           6 :         call sortAscending(ndata,WorkSpace2,WorkSpace1,Err)
     116           6 :         if (Err%occurred) then
     117             :         ! LCOV_EXCL_START
     118             :             Err%msg = PROCEDURE_NAME//Err%msg
     119             :             return
     120             :         end if
     121             :         ! LCOV_EXCL_STOP
     122           6 :         call crank(ndata,WorkSpace2,sg)
     123             : 
     124         306 :         WorkSpace1 = WorkSpace1 - WorkSpace2
     125         306 :         dStarStar = dot_product(WorkSpace1,WorkSpace1)
     126           6 :         en = ndata
     127           6 :         en3n = en**3-en
     128           6 :         aved = en3n/6.0_RK-(sf+sg)/12.0_RK
     129           6 :         fac = (1.0_RK-sf/en3n)*(1.0_RK-sg/en3n)
     130           6 :         vard = ((en-1.0_RK)*en**2*(en+1.0_RK)**2/36.0_RK)*fac
     131           6 :         dStarStarSignif = (dStarStar-aved)/sqrt(vard)
     132           6 :         dStarStarProb = erfc( real( abs(dStarStarSignif)/SQRT2 , kind=SPR ) )
     133           6 :         rho = (1.0_RK-(6.0_RK/en3n)*(dStarStar+(sf+sg)/12.0_RK))/sqrt(fac)
     134           6 :         fac = (1.0_RK+rho)*(1.0_RK-rho)
     135          12 :         if (fac > 0.0) then
     136           6 :             t = rho*sqrt((en-2.0_RK)/fac)
     137           6 :             df = en-2.0_RK
     138           6 :             rhoProb = getBetaCDF(0.5_RK*df,0.5_RK,df/(df+t**2))
     139             :         else
     140           0 :             rhoProb = 0.0
     141             :         end if
     142             : 
     143             :     contains
     144             : 
     145          12 :         subroutine crank(ndata,Array,s)
     146           6 :             use Constants_mod, only: IK, RK
     147             :             use Misc_mod, only : arth, copyArray
     148             :             implicit none
     149             :             integer(IK), intent(in) :: ndata
     150             :             real(RK), intent(out) :: s
     151             :             real(RK), intent(inout) :: Array(ndata)
     152             :             integer(IK) :: i,ndum,nties
     153          12 :             integer(IK) :: Tstart(ndata), Tend(ndata), Tie(ndata), Idx(ndata)
     154          12 :             Idx = arth(1,1,ndata)
     155         612 :             Tie = merge(1,0,Array == eoshift(Array,-1))
     156          12 :             Tie(1) = 0
     157         612 :             Array = Idx(:)
     158         612 :             if (all(Tie == 0)) then
     159          12 :                 s = 0.0
     160          12 :                 return
     161             :             end if
     162           0 :             call copyArray( pack(Idx,Tie<eoshift(Tie,1)), Tstart, nties, ndum )
     163           0 :             Tend(1:nties) = pack(Idx(:),Tie(:)>eoshift(Tie(:),1))
     164           0 :             do i = 1,nties
     165           0 :                 Array(Tstart(i):Tend(i)) = (Tstart(i)+Tend(i))/2.0_RK
     166             :             end do
     167           0 :             Tend(1:nties) = Tend(1:nties)-Tstart(1:nties)+1
     168           0 :             s = sum(Tend(1:nties)**3-Tend(1:nties))
     169          12 :         end subroutine crank
     170             : 
     171             :     end subroutine getCorrCoefSpearman
     172             : 
     173             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     174             : 
     175             : end module CorrCoef_mod ! LCOV_EXCL_LINE

ParaMonte: Plain Powerful Parallel Monte Carlo Library 
The Computational Data Science Lab
© Copyright 2012 - 2021