R/TMLG.r

In RobustAFT: Truncated Maximum Likelihood Fit and Robust Accelerated Failure Time Regression for Gaussian and Log-Weibull Case

TMLG <-
function(X,y,delta, Beta.t,sigma.t, const,tu,control) {
# TML, Gaussian case, iteratively reweighting algorithm
maxit.sigma   <- control$maxit.sigma
tol.sigma     <- control$tol.sigma
maxit.Beta    <- control$maxit.Beta
tol.Beta      <- control$tol.Beta
Maxit.TML     <- control$Maxit.TML
tol.TML.sigma <- control$tol.TML.sigma
tol.TML.Beta  <- control$tol.TML.Beta
alg.sigma     <- control$alg.sigma
nitmon        <- control$nitmon
p <- ncol(X);  Nit    <- 1; Mxt <- 10
Beta0  <- Beta.t
sigma0 <- sigma.t
mui.t  <- X %*% as.matrix(Beta.t)
rs0    <- (y-mui.t)/sigma.t
wgt    <- ww(rs0,-tu, tu)
DeltaB <- Deltas <- 100000
Nit    <- 1
while ( ( Deltas > tol.TML.sigma | DeltaB > tol.TML.Beta ) & Nit <= Maxit.TML) {
mui   <- X %*% as.matrix(Beta0)                                                            # check
if (alg.sigma ==1) {
  zs    <- TML.SigmaG(X,y,delta, sigma0, sigma.t, mui,mui.t, wgt,const,tu, tol=tol.sigma,maxit=maxit.sigma,nitmon)
  sigma <- zs$sigma; nits <- zs$nit }
if (alg.sigma !=1) {
# determine lower and upper values for sigma
  smu  <- sml  <- sigma0;  nit  <- 1 
  funu <- funl <- TML.Ave2G(X,y,delta, smu,sigma.t, mui,mui.t, wgt,tu)-const
  while(funu > 0 & nit <= Mxt) {
    smu  <- smu*1.5
    funu <-       TML.Ave2G(X,y,delta, smu,sigma.t, mui,mui.t, wgt,tu)-const
    nit  <- nit+1}
    if (nitmon) cat("smu, funu, nit",smu,funu,nit,"\n")
  nit <- 1
  while(funl < 0 & nit <= Mxt) {
    sml  <- sml*0.5
    funl <-       TML.Ave2G(X,y,delta, sml,sigma.t, mui,mui.t, wgt,tu)-const
    nit  <- nit+1}
    if (nitmon) cat("sml, funl, nit",sml,funl,nit,"\n")
# compute new sigma
  zs <- regfal(TML.Ave2G,cc=const,lower=sml,upper=smu,nint=5,tol=tol.sigma,maxit=maxit.sigma,
               X=X,y=y,delta=delta, sigma.t=sigma.t, mui=mui, mui.t=mui.t, wgt=wgt,cu=tu)
  sigma  <- zs$solution; nits <- zs$nit}
  if (nitmon) cat("sigma", nits, "\n")
# compute new Beta
  zB   <- TML.BetaG(X,y,delta,Beta0,sigma, Beta.t,sigma.t, tu, maxit=maxit.Beta,tol=tol.Beta,nitmon)
  Beta <- zB$Beta
  DeltaB <- max(abs(Beta - Beta0)); Deltas <- abs(sigma-sigma0)
  if (nitmon) cat("Nit=",Nit,"DeltaB",DeltaB,"Deltas",Deltas, "\n")
  Beta0 <- Beta; sigma0 <- sigma
  Nit   <- Nit+1}
zres <- list(Beta=Beta,sigma=sigma,Nit=Nit)
zres}