#  source("testssim.s")

#  Mar 2, 2001
#  Tests of the hypothesis in simple linear regression
#
#    model:  y = b0 + b1*x + E
#
#  Hypothesis: b1=0

#  We will consider this problem in the context of
#  M-regression.

#  The fitting will be done with rrega
#  We will consider various possibilities; that is, 
#  different weight funtions (functional forms and parameters).

#  From the set of possible tests, we will consider
#  1. a regular t test (F test)        pt
#  2. an approximate F test (FM test)  FMapp
#  3. Wald's test                      Wsimtest

#     call requirements:
#  pt        full-model residuals -- 
#            one call to rrega
#  FMapp     sighat, and residuals from both full and reduced -- 
#            two calls to rrega
#  Wsimtest  sighat, x, residuals from full-model --
#            one call to rrega

#   Outline:
#    0. initializations
#    1. loop m times
#    1.1. loop over x pattern
#    1.1.1 loop over dist of E
#    1.1.1.1 loop over b
#    1.1.1.1.1 generate data 
#    1.1.1.1.2. fit full model -- get bhat, residuals, sig
#    1.1.1.1.3. fit reduced model
#    1.1.1.1.4. call ttest, FMapp, and Wsimtest
#    1.1.1.1.5. update performance data
#    2. compute summaries

#   Initializations
library("MASS")

####### Here's a little example #######
kbend <- 1.345
alpha <-.05
m     <- 100
n     <- 10
b     <- 0
Edist <- 'dist'
xpat  <- 'pattern'
x<-c(1:n)
xconst <- rep(0,n)
nrejfm <- 0
nrejz  <- 0
nrejt  <-0
nrejza <- 0
nrejta <-0
nrejlst <- 0
nrejrlmt <- 0
nrejrlmfm <- 0

for (i in 1:m) {
y<-1+b*x+rnorm(n)

tlm<-lm(y~x)
trlm<-rlm(y~x,scale.est="proposal 2")
rlmtpvalue  <- 2*pt(-abs(summary(trlm)$coef[2,3]),n-2)
sigrlm <- summary(trlm)$sigma
trlmred  <- rrega(xconst,y, 
                method=function(u) wt.huber(u,c=kbend),fix.scale=sigrlm)
fmrlm  <- FMapp(summary(trlm)$residuals, 
      trlmred$residuals, sigrlm, 1, 0)

rrful  <- rrega(x,y,
                method=function(u) wt.huber(u,c=kbend))
bhat   <- rrful$coef[2]
resful <- rrful$residuals
sig    <- 1.483*median(abs(resful))
rrred  <- rrega(xconst,y, 
                method=function(u) wt.huber(u,c=kbend),fix.scale=sig)
resred <- rrred$residuals
fm     <- FMapp(resful, resred, sig, 1, 0)
w      <- Wsimtest(x, bhat, resful, sig, k=kbend)

if (fm$pvalue < alpha) nrejfm <- nrejfm + 1
if (w$zpvalue < alpha)nrejz  <- nrejz + 1
if (w$tpvalue < alpha)nrejt  <- nrejt + 1
if (w$zpvaluea < alpha)nrejza <- nrejza + 1
if (w$tpvaluea < alpha)nrejta <- nrejta + 1
if (summary(tlm)$coef[2,4] < alpha) nrejlst <- nrejlst + 1 
if (rlmtpvalue  < alpha)  nrejrlmt <- nrejrlmt +1
if (fmrlm$pvalue < alpha) nrejrlmfm <- nrejrlmfm +1


}




nrejfm/m
nrejz/m
nrejt/m
nrejza /m
nrejta/m
nrejlst/m
nrejrlmt/m
nrejrlmfm/m