Finding all zeros of a univariate function on the real line

Specific Domains Numerics
1

I am trying to find all zeros of

f(z) = C_m + D_m H(C_f + D_f H(z)) - z

where C_i, D_i \in \mathbb{R} (can be any sign) and H is the CDF of the standard normal. If someone wants to experiment, this can be coded in Julia as

using StatsFuns
makef(Cm, Dm, Cf, Df) = z -> Cm + Dm * normcdf(Cf + Df * normcdf(z)) - z

Now if D_m D_f \le 2 \pi one can prove uniqueness. Existence is guaranteed because f goes to \pm \infty at the edges.

What I am interested in is finding all zeros, for generic/arbitrary arguments.

When that does not hold, eg

args = (23.449772654974293, -22.97316619909447, 0.6575825044958663, -0.6582958874152536)
makef(args...)

f

Is there a robust algorithm for this? My first thought was IntervalRootFinding.jl, but I ran into an issue [that MWE has p = normcdf(z), which as a host of numerical problems on its own, see the values above).

2

ApproxFun?

2 Likes
3

Not a solution, but you should be able to restrict your search space by noting that if f = 0, then

z = C_m + D_mH

and since H \in [0,1] (since it’s a cdf), this means that the root z must be in the interval [C_m, C_m + D_m] (if D_m is positive) or [C_m + D_m, C_m] (if D_m is negative).

1 Like
4

FWIW, you could use AccessibleModels.jl to explore the problem graphically:

using AccessibleModels, IntervalSets
using GLMakie
using StatsFuns

struct makef{A, B, C, D}
    Cm::A
    Dm::B
    Cf::C
    Df::D
end

(m::makef)(z) = m.Cm + m.Dm * normcdf(m.Cf + m.Df * normcdf(z)) - z

amodel = AccessibleModel(makef(1., 1., 1., 1.),
(   (@o _.Cm) => -100..100,
    (@o _.Dm) => -100..100,
    (@o _.Cf) => -100..100,
    (@o _.Df) => -100..100,
))

fig = Figure()
obj, = SliderGrid(fig[1,1], amodel)

lines(fig[2,1], -100..100, @lift x -> $obj(x))

And if someone can add input text boxes to the example, it would be awesome.

1 Like
5

This is a great idea! I can take this further and define

y = \frac{z - C_m}{D_m}

and then my problem is

H(C_f + D_f H(C_m + D_m y)) = y

where we know that y \in [0,1] so there is a bracket already. Brent’s method from Roots.jl is pretty good on this, not as fast as Newton’s, but much more robust.

I am only struggling with precision issues where I am around the edges (y \approx 0 or y \approx 1).

Thanks, that is really helpful. Now I convinced myself that there is either a unique root or 3 roots. Life would be so much simpler if I got IntervalRootFinding.jl working on this problem.

1 Like
6

Since the interval for y is bounded you can use isolate_roots form GitHub - Joel-Dahne/ArbExtras.jl: Julia package with extra tools for Arblib.jl which will either prove to you that all roots are isolated or indicate which ones were not.

You might define normcdf as:

julia> function normcdf!(out::Arb, x::Arb)
          t = Arb(2, prec = precision(x))
          sqrt2 = Arblib.sqrt!(t, t)
          Arblib.inv!(t,t)
          Arblib.mul!(t, t, x)
          Arblib.hypgeom_erf!(out, t)
          Arblib.add!(out, out, 1)
          Arblib.div!(out, out, 2)
          out
       end;

julia> StatsFuns.normcdf(a::Arb) = (x = zero(a); normcdf!(x, a); x);

julia> normcdf(2) - normcdf(Arb(2))
[-1.3849763108389696060605952972944153160163339998695237096472e-18 +/- 3.28e-77]

note that to use isolate_roots you’d need to define normcdf also for ArbSeries, which I’m not familiar with; If you’re interested ask @joeldahne about this!