# src/ratfor/hoeffd.r In harrelfe/Hmisc: Harrell Miscellaneous

```#  ratfor -o ../hoeffd.f hoeffd.r
#
SUBROUTINE hoeffd(xx, n, p, dmat, aadmat, madmat, npair, x, y, rx, ry, rj)
IMPLICIT DOUBLE PRECISION (a-h,o-z)
INTEGER p, npair(p,p)

DO i=1, p {
np=0
DO k=1, n {
if(xx(k,i) < 1d49) np = np + 1
}
npair(i,i) = np

DO j=(i+1),p {
m = 0
DO k=1,n {
xk = xx(k,i)
yk = xx(k,j)
if(xk < 1d49 & yk < 1d49) {
m = m + 1
x(m) = xk
y(m) = yk
}
}
npair(i,j) = m
if(m > 4) {
dmat(i,j) = d
}
else dmat(i,j) = 1d50
}
}
DO i=1,p {
dmat(i,i) = 1d0/30d0
DO j=(i+1),p {
dmat(j,i) = dmat(i,j)
npair(j,i) = npair(i,j)
}
}
RETURN
END

IMPLICIT DOUBLE PRECISION (a-h,o-z)
DOUBLE PRECISION x(n), y(n), rx(n), ry(n), rj(n), maxad
#	INTEGER iwork(1)
#	CALL rank(n, x, work, iwork, rx)
#	CALL rank(n, y, work, iwork, ry)
CALL jrank(x, y, n, rx, ry, rj)
q = 0d0
r = 0d0
s = 0d0
z = n
DO i=1,n {
rxi = rx(i)
ryi = ry(i)
rji = rj(i)
q = q + (rxi-1d0)*(rxi-2d0)*(ryi-1d0)*(ryi-2d0)
r = r + (rxi-2d0)*(ryi-2d0)*(rji-1d0)
s = s + (rji-1d0)*(rji-2d0)
}
d = (q-2d0*(z-2d0)*r+(z-2d0)*(z-3d0)*s)/z/(z-1d0)/(z-2d0)/(z-3d0)/(z-4d0)
RETURN
END

#	Use C version of this which is much faster (since uses a sort)
#	SUBROUTINE rank(x, n, r)   # simple rank with midranks for ties
#	REAL*4 x(1), r(1)
#	DO i=1,n {
#	  xi=x(i)
#	  ir=2   # will be 2*rank(x(i))
#	  DO j=1,n {
#	    if(i.ne.j) {
#	      if(x(j)<xi) ir=ir+2
#	      if(x(j)==xi) ir=ir+1
#	    }
#	  }
#	  r(i)=float(ir)/2.0
#	}
#	RETURN
#	END

SUBROUTINE jrank(x, y, n, rx, ry, r) # joint rank with midranks for ties
INTEGER n
DOUBLE PRECISION x(n), y(n), rx(n), ry(n), r(n), cx, cy, ri, rix, riy, xi, yi
DO i=1,n {
xi  = x(i)
yi  = y(i)
ri  = 1d0
rix = 1d0
riy = 1d0
DO j=1,n {
if(i .ne. j) {
cx = 0d0
if(x(j) < xi) cx = 1d0
if(x(j) == xi) cx = .5d0
cy = 0d0
if(y(j) < yi) cy = 1d0
if(y(j) == yi) cy = .5d0
rix = rix + cx
riy = riy + cy
ri  = ri + cx*cy
}
}
rx(i) = rix
ry(i) = riy
r(i) = ri
}
RETURN
END
```
harrelfe/Hmisc documentation built on May 19, 2024, 4:13 a.m.