[SciPy-dev] Possible Error in Kendall's Tau (scipy.stats.stats.kendalltau)
Sturla Molden
sturla at molden.no
Wed Mar 18 18:00:28 EDT 2009
> I was thinking about how to best use the contigency-table for tau-b. I
> don't really see how it can be done without some loops. It may be easier
> to do this in Fortran.
f2py'ing something like this should work (not thoroughly tested though)...
subroutine taub(C, D, table, tau)
implicit none
intrinsic :: sqrt, sum
integer*4, intent(in) :: C, D
integer*4, dimension(C,D), intent(in) :: table
real*8, intent(out) :: tau
integer*4 :: i, j, tmp, P, Q, Tx, Ty
P = 0
Q = 0
Tx = 0
Ty = 0
!$omp parallel do &
!$omp& private(i,j,tmp) &
!$omp& reduction(+:P,Q,Tx,Ty) &
!$omp& shared(table,C,D)
do i = 1,D
do j = 1,C
tmp = table(j,i)
! count concordant pairs
if ((i .lt. D) .and. (j .gt. 1)) then
P = P + sum(table(1:j-1,i+1:D)*tmp)
end if
! count disconcordant pairs
if ((i .gt. 1) .and. (j .lt. C)) then
Q = Q + sum(table(j+1:C,1:i-1)*tmp)
end if
! count pairs tied in y
if (i .lt. D) then
Ty = Ty + sum(table(j,i+1:D)*tmp)
end if
! count pairs tied in x
if (j .lt. C) then
Tx = Tx + sum(table(j+1:C,i)*tmp)
end if
end do
end do
!$omp end parallel do
tau = real(P-Q)/sqrt(real((P+Q+Tx)*(P+Q+Ty)))
end subroutine
More information about the SciPy-Dev
mailing list