Probability that $b^2 - 4ac \geq 0$ where $a,b,c$ are normally distributed (numerical integration)

Question

I would like to determine the probability that a random quadratic polynomial has positive discriminant, where the 3 coefficients $a, b, c$ are normally distributed and independent:

That is, given $a,b,c \sim \operatorname{N}\left(0,1\right)$, what is ( a numerical approximation to 5 digits of ) $\operatorname{Pr}\left(\,{b^2 - 4 a c \geq 0}\,\right)$ ?

---

Thoughts:

We have

$$ \mathrm{Pr}\left(\,{b^2 - 4 a c \geq 0}\,\right) = \dfrac{1}{\left(\sqrt{\pi}\right)^3} \iiint_{\mathbb{R}^3} \mathbb{\large 1}_{b^{2}\ -\ 4ac\ \geq\ 0}\quad{\rm e}^{-a^2} {\rm e}^{-b^2}{\rm e}^{-c^2}\,\mathrm{d}a \,\mathrm{d}b \,\mathrm{d}c $$

This integral probably cannot be expressed explicitly, but even a numerical approximation is not so easy. Here is what I tried with SAGE, without success:

var('a,b,c') 
RR = RealField(100)
 
I = integrate(integrate(integrate(exp(-a^2), a, 0, b^2/(4*c)) * exp(-c^2),    c, 0, oo) * exp(-b^2),    b, 0, oo) print(RR(I)) #error...

Expermenting with SAGE seems to give a probability between 0.64 and 0.65 (compare the values given below if $a,b,c$ are uniformly distributed):

N = 0
T = 10^5
for _ in range(T):
    a = gauss(0, 1)
    b = gauss(0, 1)
    c = gauss(0, 1)
    if b^2 - 4*a*c >= 0:
        N += 1

print(float(N/T))

---

Other comments:

1. The idea to consider the normal distribution is that the vector $(a, b, c) / \sqrt{a^2 + b^2 + c^2}$ is uniformly distributed on the sphere, because the joint Gaussian distribution is spherically symmetric, and the property $b^2 - 4 a c \geq 0$ is invariant under rescaling.

2. Note that when $a,b,c \sim \mathrm{Unif}(-N, N)$ for some $N > 0$, the probability $\mathrm{Pr}(b^2 - 4 a c \geq 0)$ can be computed explicitly [1]: it is $$(41 + \log(64))/72 \approx 0.627207.$$ See also Probability that a quadratic polynomial with random coefficients has real roots for the case $a,b,c \sim \mathrm{Unif}(0, 1)$.

Answer 1

Suppose we knew the expected number $\Bbb E$ of real roots of a random quadratic polynomial $aX^2 + bX + c$, where all $3$ coefficients have standard normal distribution, $\mathcal{N}(0, 1)$. Then, if $\Bbb P(\Delta \geq 0)$ is the probability that the discriminant $\Delta = b^2 - 4 a c$ of such a polynomial is positive, by definition $\Bbb E = 2 \cdot \Bbb P(\Delta \geq 0) + 0 \cdot \Bbb P(\Delta < 0) = 2 \Bbb P(\Delta \geq 0)$, so $$\Bbb P(\Delta \geq 0) = \frac12 \Bbb E .$$

On the other hand, a clever argument of Edelman and Kostman that identifies a polynomial $c + b t + a t^2$ with the vector $\langle c, b, a\rangle \in \Bbb R^3$ and considers the projection of the so-called moment curve $\gamma(t) = \langle 1, t, t^2 \rangle$ in $\Bbb R^3$ to the unit sphere shows that the expected number $\Bbb E$ of real zeros is $\frac1\pi \vert\gamma\vert$, where $|\gamma|$ denotes the arc length of $\gamma$. The arc length formula gives $$ \Bbb P(\Delta \geq 0) = \frac12 \Bbb E = \frac{|\gamma|}{2 \pi} = \frac{1}{2 \pi} \int_{-\infty}^\infty \sqrt{\frac{1}{(t^2 - 1)^2} - \frac{9 t^4}{(t^6 - 1)^2}} \,dt = \frac{1}{2 \pi} \int_{-\infty}^\infty \frac{\sqrt{t^4 + 4 t^2 + 1}}{t^4 + t^2 + 1} \,dt . $$ Integrating numerically gives $$ \color{#df0000}{\boxed{\Bbb P(\Delta \geq 0) = 0.64851\ldots}} , $$ in agreement with the Monte Carlo approximation that o.p. mentions. Even with some coaxing Mathematica wouldn't compute a form for the integral in terms of special functions, though I suspect that one can express the probability in terms of elliptic integral functions.

The integral expression for $|\gamma|$ specializes a formula of Kac (1947). Both the above formula and the argument are actually given in the reference for polynomials of arbitrary degree. For cubic polynomials, the analogous computation gives the probability $$\Bbb P(\Delta \geq 0) = \frac{1}{2 \pi} \int_{-\infty}^\infty \frac{\sqrt{(t^2 + 1)^4 + 4 t^4}}{(t^4 + 1) (t^2 + 1)} \,dt - \frac12 = 0.24637\ldots .$$

Long remark When constructing a distribution on random polynomials of degree $\leq n$ in terms of their coefficients, at least for some purposes it's more natural to choose the coefficient of the degree-$k$ term to have variance $n \choose k$ times that of the leading and constant coefficients. See, for example, the discussion in $\S$3.1 of the reference.

For $n = 2$, we have $a, c \sim N(0, 1)$ and $b \sim \mathcal N(0, 2)$. For this distribution we can proceed by choosing a linear change of coordinates adapted to the geometry of the problem and a performing a simple geometric computation in those coordinates. In the coordinates $$(x, y, z) := \left(\frac{-a + c}{2}, \frac{b}{2}, \frac{a + c}{2}\right), $$

• the joint probability distribution is $$f(x, y, z) = \frac{1}{\pi^{3 / 2}} e^{-(x^2 + y^2 + z^2)}$$ and in particular is spherically symmetric, and
• the discriminant condition $\Delta \geq 0$ becomes $$x^2 + y^2 - z^2 \geq 0, $$ which is still invariant under dilation.

So, the probability that the discriminant is nonnegative is the probability that a randomly chosen point on a sphere centered at $(0, 0, 0)$ is outside the standard cone $x^2 + y^2 - z^2 = 0$, and a standard exercise gives that this probability is $$\boxed{\Bbb P(\Delta \geq 0) = \frac{1}{\sqrt{2}} = 0.70710 \ldots} .$$

The analogous computation for $n = 3$, where the coefficients of $a, b, c, d$ of $a X^3 + b X^2 + c X + d$ have variances $1, 3, 3, 1$, respectively, gives $\Bbb P(\Delta \geq 0) = \frac12 (\sqrt 3 - 1) = 0.36602\ldots$.

A. Edelman, E. Kostlan, "How many zeros of a random polynomial are real?" Bull. Amer. Math. Soc. 32 (1995), 1-37.

Answer 2

Just another presentation of the problem. Denote $$B=\frac{b+c}{\sqrt{2}}, C=\frac{b-c}{\sqrt{2}}$$ You are interested in $\Pr(X>0)$ when $X=\frac{B^2}{2}-\frac{C^2}{2}-4A^2$ where $A,B,C$ are independent $N(0,1).$ The characteristic function of $X$ is $$\varphi(t)=E(e^{itX})=\frac{1}{\sqrt{1+t^2}}\times \frac{1}{\sqrt{1-8it}}$$ with a suitable interpretation of $\sqrt{1-8it}.$ Since $\varphi$ is integrable, inversion formula applies and the density of $X$ is $$f(x)=\frac{1}{2\pi}\int_Re^{-itx}\varphi(t)dt$$ and for $T>0$ we have by Fubini $$\Pr(0<X<T)=\int_0^Tf(x)dx=\int_R\frac{1-e^{-itT}}{it}\varphi(t)dt.$$

Answer 3

It is known that the pdf $f_{ac}$of the product of two independent random variables $a$ and $c$ with standard normal distribution is (see here):

$$f_{ac}(t)=\frac{1}{\pi}K_0(|t|), \quad t\in \mathbb R$$

where $K_0(t)$ denotes the modified Bessel function of the second kind.

Then, the probability can be obtained as follows:

$$\mathbb P (b^2 \ge 4ac)=\mathbb P \big (b^2 \ge 4ac, \, ac \le 0\big)+\mathbb P \big(b^2 \ge 4ac, \, ac>0 \big)=\\ \frac{1}{2}+\mathbb P \big (b \le -\sqrt{4ac}, \, ac>0 \big)+\mathbb P \big(b \ge \sqrt{4ac}, \, ac>0\big)=\\ \frac{1}{2}+2\int_0^\infty \left (1-\Phi\left(2\sqrt{t} \right) \right )\frac{K_0(t)}{\pi}\text{d} t= \\ \frac{1}{2}+\color{blue}{\frac{2}{\pi}\int_0^\infty \left (\frac{1}{2}-\frac{1}{2}\text{erf}\left(\frac{2\sqrt{t}}{\sqrt{2}} \right) \right )K_0(t)\text{d} t}\\ \approx \color{blue}{0.6485118}$$

where $\Phi$ is the cdf of the standard normal distribution and $\text{erf}$ is the error function, with $$\Phi(x)=\frac{1}{2}+\frac{1}{2}\text{erf}\left(\frac{x}{\sqrt{2}} \right).$$ The blue integral is numerically computed here. The probability $0.6485118$ is very close to what is obtained by simulation in the OP.