[SciPy-dev] Possible Error in Kendall's Tau (scipy.stats.stats.kendalltau)

[Sturla Molden]

> 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