#  source("FMapp.s")

#  This routine computes the approximate F statistic
#  for testing 
#      B_{dred+1} = ... = B_{dful} = 0,
#  after an M-regression fit.
#  Here dred is the number of coefficients in the reduced model
#  and  dful is the number of coefficients in the full model.
#  This is equation (5.6) in Birkes and Dodge.
#  The user must do the full and the reduced fits and obtain
#  the respective residual vectors, resful and resred, prior
#  to calling FMapprox.  These fits can be obtained using
#  the S-Plus function rrega.
#
#  resful are the regular residuals from the full model fit using M
#  resred *** are residuals from the reduced model that was fit
#         using the sigma from the full model.

#  Typical usage is
#    rr  <- rreg(---[full model]---,
#                method=function(u) wt.huber(u,c=kbend))
#    sig <- 1.483*median(abs(rr$residuals))
#    rra <- rrega(---[reduced model]---, 
#                method=function(u) wt.huber(u,c=kbend),fix.scale=sig)


FMapp  <- function(resful, resred, sig, dful, dred, kbend=1.345){

n <- length(resful)
if (n != length(resred)){
   print("residual vectors not the same length")
   stop()
}
if (dred >= dful){
   print("number of terms in the full model must be greater than number in reduced model")
   stop()
}

#  Full model

mineif  <- ifelse(resful<kbend*sig,
                  resful,   kbend*sig)
fif     <- ifelse(-kbend*sig>mineif,
                  -kbend*sig, mineif)
gif     <- ifelse(0>abs(resful)-kbend*sig,
                  0,        abs(resful)-kbend*sig)
s1      <- sum(fif*fif)
strful  <- s1+2*kbend*sig*sum(gif)

# Reduced model

mineir  <- ifelse(resred<kbend*sig,
                  resred,   kbend*sig)
fir     <- ifelse(-kbend*sig>mineir,
                  -kbend*sig, mineir)
gir     <- ifelse(0>abs(resred)-kbend*sig,
                  0,        abs(resred)-kbend*sig)
strred  <- sum(fir*fir)+2*kbend*sig*sum(gir)

if (strred <= strful){
   print("sum of transformed residuals for the reduced model is not greater than STR for full model")
   stop()
}

mzero  <- sum(gif[gif==0]+1)
lambda <- n*s1/(mzero*(n-dful-1))

FM <- (strred-strful)/((dful-dred)*lambda)
pvalue <- 1-pf(FM,1,n-2)
return(FM,lambda,pvalue)
}