Finding largest root of a function in R

Question

I have a function that is data generated and not a polynomial. For most data, it only has one root, but it can also have two roots, as in in the following plot:

uniroot() tends to fail because the function has the same sign at both a small value for h and a large value for h, which I use at the interval endpoints.

Is there a way to find the greatest root of a function in R for a function that is known to go to minus infinity for large values of its argument?

As a workaround, I have placed the search interval around a large value for h, i.e. into a region on the right where I know the function has no roots, and then used the argumend extendInt="downX". I wonder however whether there is something more robust.

As requested by the comments, here is a reproducible example:

# function whose largets root h is to be found
h.root <- function(x, h) {
  n <- length(x)
  f2 <- function(y) {
    h2 <- h*h
    (sum(((x-y)^2/h2 - 1) * dnorm((x-y)/h)) / (n*h*h2))^2
  }
  f2.int <- integrate(Vectorize(f2), -Inf, Inf)$value
  (0.5 / sqrt(pi) / f2.int / n)^0.2 - h
}

# examples for problematic seeds: 45, 36
set.seed(36)
x <- rnorm(50)
h <- seq(0.05, 0.6, 0.01)
res <- h
for (i in 1:length(h)) res[i] <- h.root(x,h[i])
plot(h, res, type="l")
abline(h=0, lty=2)

print(uniroot(h.root, c(0.05, 10), x=x)$root)

Probably, rootSolve::uniroot.all would do it, but if not please come back with a workable example. — Nir Graham
– Nir Graham, Commented Jul 24, 2024 at 16:02
@nir-graham Just evaluating the function at a grid of 100 points and looking for sign changes like uniroot.all does is a possible approach, but far from robust. I have added code implementing the function and an example with random data where uniroot returns the "wrong" root. Yes, in this case uniroot.allwould work, but I am not looking for a solution that just works for normally distributed data. Maybe pracma::fzero is a way (must check the reference), but maybe someone has a solution for just the problem of finding the largest root. — cdalitz
– cdalitz, Commented Jul 24, 2024 at 16:55
If we can assume that the largest root is within 0.05 and 0.6, or other fixed bounds and that roots are at least 0.01 (or other) apart then working backwards from the upper bound to the first sign change and then running uniroot is pretty easy to do and runs fast so I would not exclude it. hseq <- seq(0.05, 0.6, 0.01); s <- h.root(x, max(hseq)); for(h in rev(hseq)) if (sign(h.root(x, h)) != sign(s)) break; uniroot(\(h) h.root(x, h), c(h, 0.6)) — G. Grothendieck
– G. Grothendieck, Commented Jul 24, 2024 at 18:39
@g-grothendieck We can assume that the root is greater than zero and less than max(x) - min(x), which may be arbitrary large, depending on the data. — cdalitz
– cdalitz, Commented Jul 24, 2024 at 20:00

ThomasIsCoding · Accepted Answer · 2024-07-26 15:03:22Z

1

Update

If you care about the largest root only, a minor change to previous code is to "locate the interval where the sign change happens", e.g.,

> idx <- max(which(rowSums(sign(embed(res, 2))) == 0))

> uniroot(h.root, c(h[idx], h[idx + 1]), x = x)$root
[1] 0.5524001

Previous Answer

Probably you can use a numerical way to solve the largest root

idx <- which(rowSums(sign(embed(res, 2))) == 0)
sols <- rowMeans(cbind(h[idx], h[idx + 1]))

and you obtain the roots (from the left to the right, where the largest is the right-most one)

> sols
[1] 0.1065 0.2545 0.5525

and you can verify the solutions

> h.root(x, sols)
[1] -7.801779e-05 -1.402256e-04 -6.189450e-05

and verify the results in the plot

plot(h, res, type = "l")
abline(h = 0, lty = 2)
points(sols, h.root(x, sols), col = "red", cex = 3)

Data

Actually you can vectorize your h.root and provide higher resolution of h (spacing by 1e-3 or even denser), e.g.,

# function whose largets root h is to be found
h.root <- Vectorize(function(x, h) {
    n <- length(x)
    f2 <- function(y) {
        h2 <- h * h
        (sum(((x - y)^2 / h2 - 1) * dnorm((x - y) / h)) / (n * h * h2))^2
    }
    f2.int <- integrate(Vectorize(f2), -Inf, Inf)$value
    (0.5 / sqrt(pi) / f2.int / n)^0.2 - h
}, vectorize.args = "h")

# examples for problematic seeds: 45, 36
set.seed(36)
x <- rnorm(50)
h <- seq(0.05, 0.6, 1e-3)
res <- h.root(x, h)

edited Jul 26, 2024 at 15:03

answered Jul 24, 2024 at 20:23

ThomasIsCoding

106k9 gold badges38 silver badges110 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

cdalitz Over a year ago

Thanks for your answer. The brute force approach has an unfavorable runtime (it is O(n*datawidth/resolution), but it can be easily speeded up by only checking for sign changes in a coarse grid and then hunting down the roots with bisection searches or with Brent's method (implemented in uniroot). AFAIK this is what uniroot.all does.

ThomasIsCoding Over a year ago

@cdalitz yes, you are right. See my update that use uniroot with refined interval to solve the largest root

Collectives™ on Stack Overflow

Finding largest root of a function in R

1 Answer 1

Update

Previous Answer

Data

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Update

Previous Answer

Data

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related