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
|
