# 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)
DOUBLE PRECISION xx(n,p), dmat(p,p), aadmat(p,p), madmat(p,p), x(n), y(n), rx(n), ry(n), rj(n), maxad
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) {
CALL hoeff(x, y, m, d, aad, maxad, rx, ry, rj)
dmat(i,j) = d
aadmat(i,j) = aad
madmat(i,j) = maxad
}
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)
aadmat(j,i) = aadmat(i,j)
madmat(j,i) = madmat(i,j)
}
}
RETURN
END
SUBROUTINE hoeff(x, y, n, d, aad, maxad, rx, ry, rj)
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
aad = 0d0
maxad = 0d0
z = n
DO i=1,n {
rxi = rx(i)
ryi = ry(i)
rji = rj(i)
ad = dabs((rji/z) - (rxi/z)*(ryi/z))
aad = aad + ad
maxad = dmax1(maxad, ad)
q = q + (rxi-1d0)*(rxi-2d0)*(ryi-1d0)*(ryi-2d0)
r = r + (rxi-2d0)*(ryi-2d0)*(rji-1d0)
s = s + (rji-1d0)*(rji-2d0)
}
aad = aad / z
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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.