Could we use uniroot() instead of optimize() for this minimization in R?

In the last line of my R code below, I use optimize() to find the df2 that minimizes the ncp_diff function.

However, I was wondering if I could uniroot() instead of optimize() for this minimization?

alpha = c(.025, .975); df1 = 3; peta = .3   # The input

f <- function(alpha, q, df1, df2, ncp){     # Notice `ncp` is the unknown
  alpha - suppressWarnings(pf(q = (peta / df1) / ((1 - peta)/df2), df1, df2, ncp, lower = FALSE))
}

ncp <- function(df2){      # Root finding: finds 2 `ncp` for a given `df2`

 b <- sapply(c(alpha[1], alpha[2]),
       function(x) uniroot(f, c(0, 1e7), alpha = x, q = peta, df1 = df1, df2 = df2)[[1]])

 b / (b + (df2 + 4))
}

ncp_diff <- function(df2, target = 0.15){
 the_ncp <- ncp(df2)
  return(abs(abs(the_ncp[2] - the_ncp[1]) - target))
 }

optimize(ncp_diff, c(0, 1000)) ## HERE can I use `uniroot()` instead of `optimize()`

1 answer

  • answered 2018-07-12 06:33 Onyambu

    alpha = c(.025, .975); df1 = 3; peta = .3   # The input
    
    f <- function(alpha, q, df1, df2, ncp){     # Notice `ncp` is the unknown
        alpha - suppressWarnings(pf(q = (peta / df1) / ((1 - peta)/df2), df1, df2, ncp, lower = FALSE))
    }
    
    ncp <- function(df2){      # Root finding: finds 2 `ncp` for a given `df2`
    
        b <- sapply(c(alpha[1], alpha[2]),
                    function(x) uniroot(f, c(0, 1e7), alpha = x, q = peta, df1 = df1, df2 = df2)[[1]])
    
        b / (b + (df2 + 4))
    }
    
    ncp_diff <- function(df2, target = 0.15){
        the_ncp <- ncp(df2)
        return((the_ncp[2] - the_ncp[1]) - target)
    }
    
    uniroot(ncp_diff, c(100, 1000)) #
    $root
    [1] 336.3956
    
    $f.root
    [1] 3.74663e-09
    
    $iter
    [1] 7
    
    $init.it
    [1] NA
    
    $estim.prec
    [1] 6.103516e-05
    

    Edit:

    In order to use the same interval of (0,1000) we can look for a way to go around a scenario where both the lower and upper values produce results on the same side of the number line. Since this is an error in r, we can go through it by tryCatch

    ncp <- function(df2){      # Root finding: finds 2 `ncp` for a given `df2`
    
      b <- sapply(c(alpha[1], alpha[2]),
                  function(x) 
        tryCatch(uniroot(f, c(0, 1e7), alpha = x, q = peta, df1 = df1, df2 = df2)[[1]],
                 error =function(e)NA ))
      if(any(is.na(b)))b= c(1,10)
    
      b / (b + (df2 + 4))
    }
    
    
    uniroot(ncp_diff, c(0, 1000)) #
    $root
    [1] 336.3956
    
    $f.root
    [1] -2.132438e-09
    
    $iter
    [1] 8
    
    $init.it
    [1] NA
    
    $estim.prec
    [1] 6.103516e-05