# SCCS @(#)frailty.controldf.s 1.3 12/02/98 # A function to calibrate the df # very empirical # Find the closest 3 points that span the target value # We know the function is monotone, so fit the function # dy = a * (dx)^p to the 3 points, where dx and dy are the distance # from the leftmost of the three points. # This method can fail near a boundary, so use step halving if things don't # go well # On input, parms$df = target degrees of freedom # parms$dfs, parms$thetas = known values (usually 0,0) # parms$guess = first guess # frailty.controldf <- function(parms, iter, old, df) { if (iter==0) { theta <- parms$guess return(list(theta=theta, done=FALSE, history=cbind(thetas=parms$thetas, dfs=parms$dfs))) } eps <- parms$eps if (length(eps)==0) eps <- .1 thetas <- c(old$history[,1], old$theta) dfs <- c(old$history[,2], df) nx <- length(thetas) if (nx==2) { #linear guess based on first two # but try extra hard to bracket the root theta <- thetas[1] + (thetas[2]-thetas[1])*(parms$df - dfs[1])/ (dfs[2] - dfs[1]) if (parms$df > df) theta <- theta * 1.5 return(list(theta=theta, done=FALSE, history=cbind(thetas=thetas, dfs=dfs), half=0)) } else{ # Now, thetas= our guesses at theta # dfs = the degrees of freedom for each guess done <- (iter>1 && (abs(dfs[nx]-parms$df) < eps)) # look for a new minimum x <- thetas y <- dfs target <- parms$df # How am I doing if ( abs( (y[nx]-target)/(y[nx-1]-target)) > .6) doing.well <- FALSE else doing.well <- TRUE ord <- order(x) if ((x[1]-x[2])*(y[1]-y[2]) >0) y <- y[ord] #monotone up else { #monotone down y <- -1* y[ord] target <- -target } x <- x[ord] if (all(y>target)) b1 <- 1 #points 1:3 are the closest then else if (all(y1 && ((target -y[b1]) < (y[b1+1] -target)))) b1 <- b1-1 } #now have the best 3 points # fit them with a power curve anchored at the leftmost one b2 <- b1 + 1:2 xx <- log(x[b2] - x[b1]) yy <- log(y[b2] - y[b1]) power <- diff(yy)/diff(xx) a <- yy[1] - power*xx[1] newx <- (log(target -y[b1]) - a)/power if (length(parms$trace) && parms$trace){ print(cbind(thetas=thetas, dfs=dfs)) cat(" new theta=" , format(x[b1] + exp(newx)), "\n\n") } list(theta=x[b1] + exp(newx), done=done, history=cbind(thetas=thetas, dfs=dfs), half=0) } }