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?
