R/unirootRobust.R

In bpcp: Beta Product Confidence Procedure for Right Censored Data

## working limit gets a limit for uniroot if the given 
## limit is +/- Inf, returns NULL if cannot find a limit
workingLimit<-function(limit,f,signTarget, start=1){
  if (abs(limit)!=Inf){
    output<- limit
  } else {
    for (i in 0:50){
      # if sign(limit)=1 keep adding
      # if sign(limit)=-1 keep substracting
      x<- ifelse(i==0,start,
                 x+ sign(limit)*2^i)
      fcurrent<- f(x)
      print(paste("x=",x," f(x)=",fcurrent))
      if (signTarget==sign(fcurrent)) break()
    }
    output<-x 
    if (signTarget!=sign(fcurrent)){
      #warning(paste0("finite limit should be more extreme than ", x))
      output<-NA
      attr(output,"limit") <- x
    } 
  }
  output
}

#f<-function(x){ -x + 2^51} 
#workingLimit(Inf,f,-1)
#workingLimit(-Inf,f,-1,start=100)



# unirootRobust still returns a result if lowerLimit=-Inf 
# or upperLimit=Inf, or if f(lower)=NA or f(upper)=NA
unirootRobust<-function(f,lower,upper){
  flo<- f(lower)
  fup<- f(upper)
  if (!is.na(flo) & !is.na(fup) ){
    if (sign(flo)==sign(fup)){
      stop("sign(f(lower))=sign(f(upper))")
    } 
    if (sign(flo)==0){
      output<- lower
    } else if (sign(fup)==0){
      output<- upper
    } else if (lower!=-Inf & upper!=Inf){
      output<- uniroot(f,c(lower,upper))$root
    } else if (lower==-Inf & upper!=Inf & fup!=0){
      # get workingLowerLimit
      signtarget<- -1*sign(fup)
      lower<-workingLimit(lower,f,signtarget)
      output<- uniroot(f,c(lower,upper))$root
      if (is.na(lower)) lower<- attr(lower,"limit")
    } else if (lower!=-Inf & upper==Inf){
      # get workingUpperLimit
      signtarget<- -1*sign(flo)
      upper<-workingLimit(upper,f,signtarget)
      output<- uniroot(f,c(lower,upper))$root
    }  else if (lower==-Inf & upper==Inf){
      # get workingLowerLimit
      signtarget<- -1*sign(fup)
      lower<-workingLimit(lower,f,signtarget)
      # get workingUpperLimit
      signtarget<- -1*sign(flo)
      upper<-workingLimit(upper,f,signtarget)
      output<- uniroot(f,c(lower,upper))$root
    }
  } else if (is.na(flo) & !is.na(fup)){
    # get workingLowerLimit
    signtarget<- -1*sign(fup)
    lower<-workingLimit(lower,f,signtarget)
    output<- uniroot(f,c(lower,upper))$root
  } else if (!is.na(flo) & is.na(fup)){
    # get workingUpperLimit
    signtarget<- -1*sign(flo)
    upper<-workingLimit(upper,f,signtarget)
    output<- uniroot(f,c(lower,upper))$root
  } else stop("both f(lower)=NA and f(upper)=NA")
  output
}


#f<-function(x){ log(x)/x + 1 }
# test 
#f(Inf)
#f(0)
#f(1)
#f(-1)
#f(-Inf)

#unirootRobust(f,lower=0,upper=1)
# unirootRobust(f,lower=0,upper=1)