A probabilistic integral $\int_{-\infty}^{\infty}e^{-x^2/2\sigma^2}\arcsin\left(1-2\left|\lfloor x\rceil-x\right|\right)\,dx$

In my probabilistic studies, a tough integral appeared. Note that $\lfloor x\rceil$ is not the floor function; it is the nearest integer function. Up to some constants, it appears in a Buffon-like situation. $$ \int_{-\infty}^{\infty}{\rm e}^{-x^2/2\sigma^{2}} \arcsin\left(1 - 2\left\vert\vphantom{\Large A}\,\left\lfloor\,x\,\right\rceil-x\,\right\vert\right) \,{\rm d}x $$

It is possible to evaluate this integral? I will be already very glad if the case $\sigma=1$ is settled. What makes this hard to me is the $\left\vert\vphantom{\Large A}\left\lfloor\,x\,\right\rceil-x\,\right\vert$ part. A natural candidate would be a Fourier expansion, but even this seems hard.

A picture of the integrand.

After this analysis, I don't think we can do much better for this integral. Thank you for your time and detailed answers.

In essence, we're interested in the expectation value of $f(x)=\arcsin\left(1-2\big|\lfloor x\rceil-x\big|\right)$ with respect to a centered Gaussian distribution:

$$\text{E}[f(x)]=\frac{1}{\sqrt{2\pi \sigma^2}}\int_{-\infty}^\infty e^{-x^2/2\sigma^2} f(x)\,dx $$

Note that $f(x)$ is a $1$-periodic function that behaves as $g(x)=\arcsin(1-2|x|)$ for $|x|\leq \frac{1}{2}$. Taking $g(x)$ to have support only on this interval, we can write $f(x)$ in a form amenable to expansion via the Poisson summation formula:

\begin{align} f(x)=\sum_{n=-\infty}^\infty g(x+n)=\sum_{k=-\infty}^\infty e^{2\pi i k x}\int_{-1/2}^{1/2} e^{-2\pi i k y}g(y)\,dy \tag{(1)} \end{align}

Note that the integral in this expression is just the $k$-th Fourier coefficient $c_k$ of $f(x)$. So we have traded the periodicity of $f(x)$ for a weighted sum of double integrals. But these all factorize into single integrals, and the integrals over $x$ are all Gaussians: \begin{align} \frac{1}{\sqrt{2\pi \sigma^2}}\int_{-\infty}^\infty \exp\left(-\frac{x^2}{2\sigma^2}+2\pi i k x\right)\,dx &=\frac{e^{-2\pi^2 k^2\sigma^2}}{\sqrt{2\pi \sigma^2}} \int_{-\infty}^\infty \exp\left(-\frac{1}{2\sigma^2}\left(x-2 \pi i k\sigma^2 \right)^2\right)\,dx\\ &=e^{-2\pi^2 k^2\sigma^2} \end{align} Note the convenient disappearance of the normalization.

Thus we may write the expectation value as $$\text{E}[f(x)] = \sum_{k=-\infty}^\infty c_k \, e^{-2\pi^2 k^2\sigma^2} $$ Noting the quadratic dependence on $k$, this series is analogous to a Jacobi theta function except weighted by the Fourier coefficients $c_k$. To proceed further we need these coefficients: \begin{align} c_k &=2\,\text{Re}\left[\int_{0}^{1/2} e^{2\pi i k y} \arcsin(1-2y)\,dy\right] &(\text{symmetry of integrand})\\ &=\text{Re}\left[e^{2\pi i k}\int_{0}^{\pi/2} e^{-2\pi i k\sin\phi} \phi \cos\phi\,d\phi\right] &(\sin\phi=1-2y) \end{align}

To proceed further, we use a Jacobi-Anger expansion to write the inner integral as $$\int_{0}^{\pi/2} e^{-2\pi i k\sin\phi} \phi \cos\phi\,d\phi =\sum\limits_{n=-\infty}^{\infty} J_n(-2\pi k) \int_{0}^{\pi/2} e^{i n \phi} \phi \cos\phi\,d\phi $$

This integral is certainly doable, but right now I don't know a nice formula for it and so halt for now. Even if is relatively nice, though, one still appears destined to end up with a doubly-infinite summation over Bessel functions and quadratic exponents. That does not seem encouraging.

UPDATE) Following the discussion in comments with @rajb245, I realized that his insight regarding the presence of the zeroth Struve H-function could be revealed directly. First, using symmetry along with the fact that $\cos(\pi k -\pi k u)=(-1)^k \cos{\pi k u}$, we rewrite $c_k$ using the substitution $y=1-2u$ to obtain $$c_k =2\int_{0}^{1/2} \cos{2\pi k y} \arcsin(1-2y)\,dy =(-1)^k \int_{0}^{1} \cos(\pi k u) \arcsin u\,du.$$ If $k=0$ then we may compute $c_0=\pi/2-1$. Otherwise, we can use the approach of this Ron Gordon answer, and integrate by parts:

\begin{align} (-1)^k \int_{0}^{1} \cos(\pi k u) \arcsin u\,du &= (-1)^k \left[\frac{1}{\pi k}\sin(\pi k u) \arcsin u-\frac{1}{\pi k}\int \frac{\sin(\pi k u)}{\sqrt{1-u^2}}\,du\right]_0^1 \end{align} The boundary term vanishes, and for the second term we recall the following integral representation of the zeroth Struve H-function: $$\mathbf{H}_0(t)=\frac{2}{\pi}\int_0^1 \frac{\sin{t u}}{\sqrt{1-u^2}}\,du$$

Therefore $c_k=\dfrac{(-1)^{k-1}}{2k}\mathbf{H}_0(2\pi k)$ for $k\neq 0$. Note that $c)k=c_{-k}$ (consistent with $g(x)$ being even) since $\mathbf{H}_0(t)$ is an odd function.

We may now return to E$[f(x)]$ and write

$$\text{E}[f(x)] =\left(\frac{\pi}{2}-1\right)-\sum_{k=1}^\infty \dfrac{(-1)^{k}}{k}\mathbf{H}_0(2\pi k)e^{-2\pi^2 k^2\sigma^2}$$

This represents my point of conclusion, since I haven't any clue whether this series can be simplified further. But this present form is quite useful: We have a constant leading term (i.e. the mean value of $\arcsin(1-2|x|)$ over $[-1/2,1/2]$) followed by a series of corrections which are each exponentially smaller than the last. Hence even taking the first correction will provide an excellent approximation. (It would probably be wise to validate this expression, either numerically or by computing the $\sigma\to0$ and $\sigma\to\infty$ limits.)

Credit goes to Semiclassical for doing most of the work on this one. My answer is only different in that I don't directly appeal to the Poisson summation formula, but just plug the Fourier series into the integral directly; a tenable sum results. Consider the integral from Semiclassical,

$$I=\frac{1}{\sqrt{2\pi \sigma^2}}\int_{-\infty}^\infty e^{-x^2/2\sigma^2} f(x)\,dx $$ with $f(x)=\arcsin\left(1-2\big|\lfloor x\rceil-x\big|\right)$. As the OP noted, Fourier expansion is a possibility since $f(x)$ is periodic with period $1$. Let us try this approach. The expansion is $$ f(x) = \sum_{k=-\infty}^\infty c_k e^{2\pi i k x}, $$ Where the coefficients are given by the integral $$ c_k = \int_{-1/2}^{1/2} f(x) e^{-i2\pi kx}\;dx. $$ We can calculate the coefficents. We'll return to this. For now, insert the Fourier expansion into the original integral $$ I=\frac{1}{\sqrt{2\pi \sigma^2}}\int_{-\infty}^\infty e^{-x^2/2\sigma^2} \sum_{k=-\infty}^\infty c_k e^{2\pi i k x}\,dx . $$ Switching integration and summation, $$ I=\frac{1}{\sqrt{2\pi \sigma^2}}\sum_{k=-\infty}^\infty c_k\int_{-\infty}^\infty e^{2\pi i k x-x^2/2\sigma^2}\,dx. $$ So long as $\sigma$ is real, that integral can be calculated for any $k$ in closed form as $$ \int_{-\infty}^\infty e^{2\pi i k x-x^2/2\sigma^2}\,dx = \sqrt{2 \pi\sigma^2 } e^{-2 \pi ^2 k^2 \sigma ^2}, $$ which leaves us with (after canceling the $\sqrt{2 \pi\sigma^2 }$) $$ I=\sum_{k=-\infty}^\infty c_k e^{-2 \pi ^2 k^2 \sigma ^2} $$

Now returning to the $c_k$ coefficients, we have $$ c_k = \int_{-1/2}^{1/2} \arcsin\left(1-2\big|x\big|\right) e^{-i2\pi kx}\;dx. $$ I'm not clever enough to exploit the symmetries here to take a shortcut, so let's just bruteforce this by splitting the integral up into two pieces, breaking the interval up at $x=0$: $$ c_k = c_k^{-} + c_k^{+} $$ with $$ c_k^{-} = \int_{-1/2}^{0} \arcsin\left(1-2\big|x\big|\right) e^{-i2\pi kx}\;dx. $$ $$ c_k^{+} = \int_{0}^{1/2} \arcsin\left(1-2\big|x\big|\right) e^{-i2\pi kx}\;dx. $$ On the negative interval, $\big|x\big|=-x$. On the positive one, $\big|x\big|=x$. The two expressions are then $$ c_k^{-} = \int_{-1/2}^{0} \arcsin\left(1+2x\right) e^{-i2\pi kx}\;dx. $$ $$ c_k^{+} = \int_{0}^{1/2} \arcsin\left(1-2x\right) e^{-i2\pi kx}\;dx. $$ On the negative interval, we'll take the substitution $x=\frac{\sin(\phi)-1}{2}$, $dx=\frac{1}{2}\cos(\phi)d\phi$, and the interval maps as $[-1/2,0]\rightarrow[0,\pi/2]$.

On the positive interval, we'll use $x=\frac{1-\sin(\phi)}{2}$, $dx=-\frac{1}{2}\cos(\phi)d\phi$, and the interval maps as $[0,1/2]\rightarrow[\pi/2,0]$. Using these substitutions now give us $$ c_k^{-} = \frac{1}{2}\int_{0}^{\pi/2} e^{i \pi k (\sin (\phi )-1)}\phi \cos (\phi )\;d\phi. $$ $$ c_k^{+} = -\frac{1}{2}\int_{\pi/2}^{0} e^{i \pi k (1-\sin (\phi ))}\phi \cos (\phi )\;d\phi. $$ Doing these integrals results in expressions involving Struve functions and Bessel functions of order zero: $$ c_k^{-} = -\frac{i \left(1-e^{-i \pi k} (J_0(k \pi )+i \pmb{H}_0(k \pi ))\right)}{4 k}. $$ $$ c_k^{+} = \frac{i \left(1-e^{i \pi k} (J_0(k \pi )-i \pmb{H}_0(k \pi ))\right)}{4 k}. $$ Now summing these two results and simplifying gives the Fourier coefficients: $$ c_k = c_k^{-}+c_k^{+} = \frac{J_0(k\pi)\sin(k\pi)-\pmb{H}_0(k\pi)\cos(k\pi)}{2 k} $$ And because $k$ is an integer, the sines and cosines simplify, leaving $$ c_k = -\frac{(-1)^k\pmb{H}_0(k\pi)}{2 k} $$ Note that the $c_0$ average value coefficient is $c_0=\frac{\pi}{2}-1$, and must be calculated by taking $\lim_{k\rightarrow 0} c_k$ since the expression has a removable singularity there.

So now, if we have a way to calculate the zero order Struve function, we have the Fourier series coefficients to plug into the summation given above for the integral $I$. Note that you can expect the summation to converge pretty fast to the ``true'' value since the terms in the sum decay like a Gaussian tail in $k$. The convergence is hastened by larger values of $\sigma$ as well, so only a few terms should give accurate results as $\sigma$ is increased.

For fun, I calculated $c_k$ for $k=0,1,2\ldots 10$ to sixteen digits, about machine precision for double precision floating point numbers on most modern computers. You get $$ 0.5707963267948966, 0.2589127103423059, 0.03248993478100901,\\ 0.04215072581916797, 0.01378829513227803, 0.01838120431720631, \\ 0.008091939971224002, 0.01070468107470534, 0.005487219243638216, \\ 0.007169482259113914, 0.004040382366534786 $$ I've plotted the corresponding finite Fourier approximation of $f(x)$, and the true $f(x)$ below, on one period.

It can be seen that with $11$ real terms, you aren't that close to the original integrand; because of the derivative blow-ups in $f(x)$, the Fourier series converges slowly to $f$. However, as we'll see below, $11$ terms is overkill for the integral with respect to a Gaussian weight.

You can show that $c_{-k}=c_k$, so you don't have to calculate the negative coefficients. The exponential term has the same symmetry, so you only have to sum over the positive values and double the result (taking care with $c_0$). The finite truncation of our integral can be written. $$ I=c_0 + 2\sum_{k=1}^N c_k e^{-2 \pi^2 k^2 \sigma ^2} $$ We can see that the exponential factor in the first term in the summation (for $\sigma=1$) is $e^{-2 \pi^2}\approx 3\times 10^{-9}$, and subsequent terms are even smaller. Essentially, I think unless $\sigma$ is on the order of $\frac{1}{\pi\sqrt{2}}$ or smaller, the entire sum is dominated by the first term, $$ I\approx c_0 = \frac{\pi}{2}-1. $$

I have numerical calculated the integral directly and verified the result. So in the limit $\sigma>>\frac{1}{2\sqrt{\pi}}$, $I=\frac{\pi}{2}-1$. At the other end, when $\sigma$ becomes very small, the exponential factor is approximately constant, and the integral would be approximately the sum of all the Fourier coefficients; I have a hunch that it sums to $\pi/2$ but haven't proven that. Alternately, you can directly use a lot of terms in the Fourier expansion for $I$. Finally, around $\sigma\approx\frac{1}{2\sqrt{\pi}}$, you'd probably not need more than 5 or 6 terms in the sum, because $\exp(-k^2)$ is order of machine precision by $k=6$.

Hmm let's see if I'm correct, I did not check it entirely but it should work this way. Maybe there are some issues at $(\star)$ like missing constants.

We know about $\lfloor x \rceil$:

$$ \forall n\in\mathbb{N}\forall x\in \left[n-\frac{1}{2},n+\frac{1}{2}\right] \quad\quad\quad \left\lfloor x\right\rceil= n \tag{$\alpha$}$$

So let's split up that integral:

$$ I:=\int\limits_{-\infty}^\infty \exp\left(-\frac{x^2}{2\sigma^2}\right)\arcsin\left(1-2\big|\lfloor x\rceil-x\big|\right) \text{d}x = {\sum_{n=-\infty}^\infty\int\limits_{n-\frac{1}{2}}^{n+\frac{1}{2}} \exp\left(-\frac{x^2}{2\sigma^2}\right)\arcsin\underbrace{\left(1-2\big|\lfloor x\rceil-x\big|\right)}_{=:A(x)} \text{d}x } $$

Now let's find a nice way of expressing $A$ throught $(\alpha)$: $$ \forall n\in\mathbb{N}\forall x\in \left[n-\frac{1}{2},n+\frac{1}{2}\right]\quad\quad\quad A(x) = 1-2\big|n-x|\le 1 $$

In fact we know even more: $$ A(x) = 1-2\big|n-x|\le 1 =1-2\big|x-n\big|=\left\{\begin{matrix} x-n+\frac{1}{2} & x\le n \\ n+\frac{1}{2}-x & x\ge n\end{matrix}\right. \tag{$\beta$}$$

So lets split the integral even further:

$$ I = {\sum_{n=-\infty}^\infty\int\limits_{n-\frac{1}{2}}^{n} \exp\left(-\frac{x^2}{2\sigma^2}\right)\arcsin\left(x-n+\frac{1}{2}\right) \text{d}x }+ {\int\limits_{n}^{n+\frac{1}{2}} \exp\left(-\frac{x^2}{2\sigma^2}\right)\arcsin\left(n+\frac{1}{2}-x\right) \text{d}x }$$

Now this is amazing because we can substitue $y=x-n+\frac{1}{2}$ and $z=n+\frac{1}{2}-x$ with $\text{d}x=\text{d}y = -\text{d}z$: $$ I= {\sum_{n=-\infty}^\infty\underbrace{\int\limits_{0}^{\frac{1}{2}} \exp\left(-\frac{\left(-y+n-\frac{1}{2}\right)^2}{2\sigma^2}\right)\arcsin\left(y\right) \text{d}y }_{B(n)}}- {\underbrace{\int\limits_{\frac{1}{2}}^{0} \exp\left(-\frac{\left(-n-\frac{1}{2}+z\right)^2}{2\sigma^2}\right)\arcsin\left(z\right) \text{d}z }_{C(n)}}$$

Now we can use the linearity of the integral that: $$ \sum_{n=-\infty}^\infty B(n) = \int\limits_{0}^{\frac{1}{2}} \arcsin\left(y\right) \sum_{n=-\infty}^\infty\exp\left(-\frac{\left(y-n+\frac{1}{2}\right)^2}{2\sigma^2}\right)\text{d}y $$

Edit: this is Wrong

It's just an inequality

Now we know for a monotonous function that: $$ \sum_{n=-\infty}^\infty B(n) =\int\limits_{0}^\frac{1}{2}\arcsin(y)\sum_{n=-\infty}^\infty\exp\left(-\frac{\left(y-n+\frac{1}{2}\right)^2}{2\sigma^2}\right) \text{d}y = {\int\limits_{0}^\frac{1}{2}\arcsin(y)\int\limits_{-\infty}^\infty \exp\left(-\frac{\left(y-n+\frac{1}{2}\right)^2}{2\sigma^2}\right) \text{d}n \text{d}y} =\sqrt{2\pi\sigma^2}\int\limits_{0}^\frac{1}{2}\arcsin(y)\text{d}y \tag{$\star$}$$

So the same is true for $C(n)$. So we get:

$$ I = 2\sqrt{2\pi\sigma^2}\int\limits_{0}^\frac{1}{2} \arcsin(x)\text{d}x $$