Compute r(nu_1,nu_2) from Marden and Perlman (1982)

marden_perlman_1982_A_1.R

# These functions are a translation into R of Marden and 
# Perlman's (1982) FORTRAN code in Appendix A1 (p. 178). 
# Marden, J. I., & Perlman, M. D. (1982). THE MINIMAL COMPLETE CLASS OF 
# PROCEDURES FOR COMBINING INDEPENDENT NONCENTRAL F-TESTS. 
# In S. S. Gupta & J. O. Berger (Eds.), Statistical Decision Theory and 
# Related Topics III (pp. 139–181). Academic Press. 
# https://doi.org/10.1016/B978-0-12-307502-4.50014-3

rstar = function(nu, mu, eps = 5e-10, maxiter = 25){
  minrstar = max(.5, nu/(nu+mu))
  y0 = nu/2 - .5
  y1 = y0 + 1
  d0 = DT(nu,mu,y0)
  d1 = DT(nu,mu,y1)
  for(k in 1:maxiter){
    yN = y0 - d0 * (y1-y0) / (d1-d0)
    dN = DT(nu,mu,yN)
    if(is.nan(dN)) return(NA_real_)
    if(abs(dN) < eps) break
    if(k == maxiter) return(NA_real_)
    y1=y0
    y0=yN
    d1=d0
    d0=dN
  }
  t0 = T0(nu,mu,yN)
  if(t0>1 | t0<minrstar) return(NA_real_)
  t0
}

rstar_v = Vectorize(rstar,c("nu","mu"))

zeroF0 = function(m, n, y){
  zF0 = 1
  for(k in 1:m){
    zF0 = zF0 * (-y) * (-k)/( (m-k+1)*(n+m-k) ) + 1
  }
  return(zF0)
}

T0 = function(nu, mu, y){
  w = nu/2
  z = w + mu/2
  f0 = zeroF0(mu/2,nu/2,y)
  f1 = zeroF0(mu/2,nu/2+1,y)
  f2 = zeroF0(mu/2,nu/2+2,y)
  1 - y * ((z/w)*f1/f0-(z + 1) * f2/((w + 1)*f1))
}

DT = function(nu, mu, y){
  w = nu/2
  z = w + mu/2
  f0 = zeroF0(mu/2,nu/2,y)
  f1 = zeroF0(mu/2,nu/2+1,y)
  f2 = zeroF0(mu/2,nu/2+2,y)
  f3 = zeroF0(mu/2,nu/2+3,y)
  dT = -( (z/w) * f1/f0 - (z + 1) * f2/( (w + 1)*f1) )
  a = z*(z +1)*f2/(w*(w+1)* f0)
  a = a - z^2 * f1^2 / ( w^2*f0^2 )
  a = a - (z + 1) * (z + 2) * f3 / ((w + 1)*(w + 2)*f1)
  a = a + (z + 1) * (z + 1) * f2^2 / ((w + 1)*(w + 1) * f1^2)
  dT - y*a
}

# Works
rstar(2,10)
# Fails; according to them, this should be .8125
rstar(1,1)  

# Reproduce their table (works)
n = c(seq(2,20,2), seq(40,100,20))
outer(n,n,rstar_v)