Numerical Integration Problems I want to analyze the following integral changing $\sigma$ over the interval $[0,\infty)$
$$\frac{1}{\sigma}\int_{0}^{1}\frac{(y-1/2)^2}{y(1-y)}\exp\bigg\{-\frac{\big[\log\big(\frac{y}{1-y}\big)\big]^2}{2\sigma^2}\bigg\}dy$$
For large $\sigma$, the adaptive quadrature method has numerical problems (probably due to the singularity of the integrand at $y=1$). Is there a way to evaluate it?
I've tried also to study its equivalent form with the same method, namely
$$\frac{1}{\sigma}\int_{0}^{\infty}\frac{(t-1)^2}{4t(t+1)^2}\exp\bigg\{-\frac{(\log t)^2}{2\sigma^2}\bigg\}dt$$
where $t=\frac{y}{1-y}$ but I still get numerical problems and find huge differences between the two integrals even for small values of $\sigma$. 
 A: I am not sure whether you're testing a numerical routine in particular, but this has an analytic solution. Call
$$
I(\sigma) = \frac{1}{4\sigma}\int_0^{+\infty}{\rm d}t~ \frac{(t-1)^2}{t}\exp\left[ -\frac{\ln^2 t}{2\sigma^2}\right] \tag{1}
$$
Before we continue, note that
\begin{eqnarray}
\int_{-\infty}^{+\infty}{\rm d}x~ e^{\alpha x}e^{-x^2/2\sigma^2} &=& e^{-\alpha^2\sigma^2/2}\int_{-\infty}^{+\infty}{\rm d}x~ e^{-(x -2\alpha\sigma^2)^2/2\sigma^2} \\
&=& \sqrt{2\pi \sigma^2}e^{-\alpha^2\sigma^2/2} \tag{2}
\end{eqnarray}
Now, going back to the original problem, call $u = \ln t$, therefore
\begin{eqnarray}
I(\sigma) &=& \frac{1}{4\sigma}\int_{-\infty}^{+\infty}{\rm d}u~ (e^u - 1)^2e^{-u^2/2\sigma^2}\\
&=&\frac{1}{4\sigma}\int_{-\infty}^{+\infty}{\rm d}u~ (1 -2 e^u + e^{2u})e^{-u^2/2\sigma^2} \\
&\stackrel{(2)}{=}& \frac{\sqrt{2\pi}}{4}\left[1 - e^{\sigma^2/2} + e^{2\sigma^2} \right] \tag{3}
\end{eqnarray}
A: For the first integral, you will want to use tanh-sinh quadrature. The quadrature is designed to handle singularities at the endpoints, which cannot in general be handled well by low-order Gauss-Konrod.
For the second integral, split the integral into $[0, 1]$ and $[1, \infty)$, use tanh-sinh quadrature on the first and exp-sinh on the second.
An implementation of tanh-sinh, exp-sinh, and sinh-sinh quadrature in arbitrary precision is given here.
A: Just to expand on the power of tanh-sinh, I've tried it with the tanh_sinh Python package.
It's important to realize that the integration points are very close to the boundary, so to integrate a finite domain, you'd integrate the first half, flip it over, and integrate the second. Since our integrand is symmetric around 0.5, it's sufficient to integrate to there and double the result:
import numpy as np
import tanh_sinh


sigma = 2.0


def f(y):
    return (
        (y - 0.5) ** 2
        / y
        / (1 - y)
        * np.exp(-(np.log(y / (1 - y)) ** 2 / 2 / sigma ** 2))
        / sigma
    )


val, err = tanh_sinh.integrate(f, 0.0, 0.5, 1.0e-10)
print(2 * val)

0.24708742951175125

And indeed, this is the correct result according to WolframAlpha.
