Possible ways to speed up slow integration?

asked 2023-07-10 23:46:45 +0200

Stew gravatar image

I am working on a problem which requires calculating the following integral: $\int_{-\infty}^{\infty}\phi(x)\Phi(x+h)^{n-1}dx$ where $\phi(x)$ and $\Phi(x)$ are the pdf and cdf of the standard Gaussian respectively. Note that $h$ and $n$ are variables where $h > 0$, $n > 0$, and $n$ is an integer.

In particular, I have the following Sage code:

x, h, n = var('x, h, n')
assume(n > 0)
assume(n, 'integer')

def pdf(x):
    return exp(-x^2/2)/sqrt(2*pi)
def cdf(x):
    return (1/2)*(1 + erf(x/sqrt(2))
f = pdf(x) * (cdf(x + h)^n)

algo = 'maxima'
result = f.integrate(x, -oo, oo, algorithm=algo)
show(result)

When I run this code, Sage just hangs and does nothing for a long time. I have even tried other algorithms such as mathematica_free and giac, but the former returned "no outputs found in the answer from Wolfram Alpha" while the latter returned "Limit: Max order reached or unable to make series expansion". I suspect that this might be a hard function to integration, but it could also be that I made an error in the code above.

Is there an error? If not, how might I get Sage to speed up the integral computations?

edit retag flag offensive close merge delete

Comments

I checked that neither Fricas, Sympy, the Wofram engine nor the Rubi integrator (used througt the Wolfram engine) can do anything more than moving some multoiplicative constant out of the integrand. Giac "never" returns (meaning that I lost patience after about 10 minutes...).

You might as well give an name to this function and derive its properties relevant to your problem in order to create a symbolic function. That would be Yet Another special function....

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2023-07-11 22:09:10 +0200 )edit

Thanks for the feedback! I actually saw this integral in Owen's A Table of Normal Integrals in which he claims that it is equivalent to some multi-variate normal cdf. I just wanted to see if this was indeed the case. Could you elaborate a bit more on the process of creating a symbolic function?

Stew gravatar imageStew ( 2023-07-12 00:53:50 +0200 )edit

Could you elaborate a bit more on the process of creating a symbolic function?

I could paraphrase function? (and the related sage.symbolic.function_factory.function), but you'll probably benefit more from your reading.

Similarly, reading the source code vor some symbolic function will show you what can and cannot be done.

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2023-07-13 08:10:52 +0200 )edit