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

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.74663e09 $iter [1] 7 $init.it [1] NA $estim.prec [1] 6.103516e05
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.132438e09 $iter [1] 8 $init.it [1] NA $estim.prec [1] 6.103516e05