# Double integral involving the normal CDF

I need to compute (or best approximate?) the following integral

$$\int_0^\infty \int_0^\infty (1 + \alpha u)^{-1}(1 + v)^{-1} \Phi\left(\frac{\beta}{\sqrt{\gamma + uv}}\right) \text{d}u \text{d}v,\qquad \gamma > 0, \quad \beta \in \mathbb{R}, \quad \alpha \in \mathbb{N}^*,$$ where $$\Phi(\cdot)$$ is the standard normal cumulative distribution function and $$\alpha$$ is typically a large natural number. What strategy would you advise for this?

I don't know if this is helpful but I can also reformulate my problem (up to multiplicative constants) as computing the following expectation:

$$\text{E}\left\{ \Phi\left(\frac{\beta}{\sqrt{\gamma + X^2Y^2}}\right)\right\},\qquad \gamma > 0, \quad \beta \in \mathbb{R},$$ where $$X \sim \text{Half-Cauchy}(0,1), \qquad Y \sim \text{Half-Cauchy}(0,\alpha^{-1/2}), \quad \alpha \in \mathbb{N}^* \, (\text{large}).$$

EDIT: I made a mistake in the first integral, it should be:

$$\int_0^\infty \int_0^\infty (1 + \alpha x^2)^{-1}(1 + y^2)^{-1} \Phi\left(\frac{\beta}{\sqrt{\gamma + x^2y^2}}\right) \text{d}x \text{d}y,\quad \gamma > 0, \beta \in \mathbb{R}, \alpha \in \mathbb{N}^*\, \text{(large)},$$ which now corresponds to the above expectation.

EDIT 2: Another way to solve my problem would be to characterise the difference (hopefully small?) $$\text{E}\left\{ \Phi\left(\frac{\beta}{\sqrt{\gamma + X^2Y^2}}\right)\right\} - \Phi\left(\frac{\beta}{\sqrt{\gamma}}\right),\qquad \gamma > 0, \quad \beta \in \mathbb{R},$$ would it be sensible to approximate the difference $$\Phi\left(\frac{\beta}{\sqrt{\gamma + X^2Y^2}}\right) - \Phi\left(\frac{\beta}{\sqrt{\gamma}}\right),$$ with a Taylor series in zero or is this unreasonable given that the half-Cauchy distributions are heavy-tailed?

• I would advise studying the asymptotic behavior, because this integral diverges.
– whuber
Mar 8, 2019 at 12:13
• Thanks for your point, I actually made a mistake in the first integral, very sorry about this. I have edited the question. Mar 8, 2019 at 12:48
• I think it would help to present the original version of this problem. This integral does not seem straightforward to evaluate using any other approach than brute force numerical quadrature, and I suspect that you have made a reformulation which is not necessarily helpful. Mar 8, 2019 at 16:19
• You are probably right, in fact the original problem corresponds to that of the expectation above, with X and Y having their respective half-Cauchy distribution... I should probably have stated the problem directly in this way. Does this make you think of any other strategy? Mar 9, 2019 at 6:46
• A more promising approach is to integrate by parts in order to introduce a Gaussian into the integrand. How to proceed from there depends on the size of $\beta/\sqrt{\gamma}.$
– whuber
Mar 9, 2019 at 15:43

Under the conditions that $$\alpha, \gamma, u,$$ and $$v$$ are all positive,

$$\beta/\sqrt{\gamma + uv} \ge \min(\beta/\sqrt{\gamma}, 0) = \delta \gt -\infty.$$

Therefore, because $$\Phi$$ is a CDF for a distribution supported on $$(-\infty,\infty),$$ $$\Phi\left(\frac{\beta}{\sqrt{\gamma + uv}}\right)\ge \Phi(\delta) = \epsilon \gt 0.$$

Consequently

\eqalign{ &\int_0^\infty \int_0^\infty (1+\alpha u)^{-1}(1+v)^{-1} \Phi\left(\frac{\beta}{\sqrt{\gamma + uv}}\right) \mathrm{d}u\mathrm{d}v \\ &\ge \epsilon\int_0^\infty \int_0^\infty (1+\alpha u)^{-1}(1+v)^{-1} \mathrm{d}u\mathrm{d}v \\ &= \lim_{M\to\infty}\lim_{N\to\infty}\epsilon\int_0^M (1+\alpha u)^{-1}\mathrm{d}u\int_0^N (1+v)^{-1} \mathrm{d}v \\ &=\frac{\epsilon}{\alpha} \lim_{M\to\infty}\lim_{N\to\infty} \log(1 + M\alpha)\log(1 + N), }

which diverges to $$+\infty.$$

• Thanks, this is a very clear argument. For the "edited" integral, it seems that the argument doesn't apply anymore. Mar 8, 2019 at 12:50
• I will let this answer stand because (a) the changes you made to the question translate in an obvious way to comparable changes in this analysis and (b) it sheds useful light on issues of numerical integration, which almost surely is going to be the solution. In particular, the convergence of each integral is only $O(1/M)$ and $O(1/N),$ which suggests preliminary manipulations to concentrate the mass of the integral into a smaller region would be useful if you want efficiency or accuracy.
– whuber
Mar 8, 2019 at 16:31
• Thanks @whuber, that's helpful, I will try this. In the meantime, I added a second edit, where I give another formulation of my problem. Mar 9, 2019 at 7:13

If approximation is OK, you could simulate the expected value through an average. A law of large numbers implies that the average converges to the expected value - provided the latter exists, which is however not obvious to me. I did simulate a few runs, though, and got very similar results each time, which is not indicative of an issue with heavy tails.

library(LaplacesDemon)

draws <- 1e6
alpha <- 3
gamma <- 1
beta <- 1
x <- rhalfcauchy(draws, scale=1)
y <- rhalfcauchy(draws, scale=1/sqrt(alpha))

mean(pnorm(beta/sqrt(gamma+x^2*y^2)))

• Thanks a lot for your reply. Yes I could do this, but in fact I would need an analytic approximation not a numeric one (in particular, alpha, beta and gamma are generic)... Mar 8, 2019 at 10:26
• OK, I hope somebody else could weigh in. Of course, you could run the code for any values of the parameters you are interested in. Mar 8, 2019 at 10:27
• Yes, sure, the problem is that I really need something analytical unfortunately. But thanks very much for your help! Mar 8, 2019 at 10:32