Incorrect evaluation of integral involving a DiracDelta, whose argument has infinitely many zeros

Question

Bug introduced in 9.0 or earlier and persisting through 12.0; reported as CASE:2323185

---

Let's say $X>0$ is a random variable with probability density $p_X(x)={\rm e}^{-x}$. Define the random variable $Y=\sin(X)$. From the transformation theorem for probabilities we know that its probability density can be obtained as the integral \begin{equation} p_Y(y)=\int_0^\infty{\rm d}x\;p_X(x)\;\delta\big(y-\sin(x)\big) = \int_0^\infty{\rm d}x\;{\rm e}^{-x}\;\delta\big(y-\sin(x)\big) \ . \end{equation} I was disappointed to see that Mathematica (8.0.1.0) misbehaves badly with this integral. The line

Integrate[Exp[-x]*DiracDelta[y - Sin[x]], {x, 0, Infinity},
          Assumptions -> {y > 0, y < 1}]    

is evaluated to $\exp\{-\arcsin(y)\}/\sqrt{1-y^2}$, an answer one gets from the rather naive substitution $z=\sin(x)$, ${\rm d}z=\cos(x)\,{\rm d}x=\sqrt{1-z^2}\,{\rm d}x$, combined with a daring interpretation of $\sin(\infty)$ in the upper integral boundary. I haven't yet had the patience to work out what the correct answer is, but it is relatively easy to see that $p_Y(0^+)=(1-{\rm e}^{-\pi})^{-1}\approx 1.045$, whereas Mathematica's answer would give $p_Y(0^+)=1$. One can also just sample the distribution and plot its histogram, which will show that Mathematica's answer$-$while looking qualitatively similar$-$is quantitatively off.

Am I using Mathematica incorrectly here? Am I expecting too much? Clearly, the trouble is that the argument of the $\delta$-function has infinitely many zeros, and so the answer will involve an infinite sum. Would I have to tell this Mathematica explicitly? Or does it really blunder ahead and not worry about these things? If the latter is the case, I'm getting quite nervous about using the $\delta$-function at all$-$except for simple cases where its argument can be uniquely inverted.

Any thoughts?

Answer 1

To elaborate on the Comment by @Jens, Integrate probably is mishandling this integrand, because ArcSin has branch cuts. The correct solution can be obtained by observing

int1 = Integrate[Exp[-x]*DiracDelta[y - Sin[x]], {x, 0, Pi/2}, Assumptions -> {0 < y < 1}]
(* 1/(E^ArcSin[y]*Sqrt[1 - y^2]) *)

int2 = Integrate[Exp[-x]*DiracDelta[y - Sin[x]], {x, Pi/2, Pi}, Assumptions -> {0 < y < 1}]
(* E^(-Pi + ArcSin[y])/Sqrt[1 - y^2] *)

Corresponding integrals in the ranges {x, Pi, 3 Pi/2} and {x, 3 Pi/2, 2 Pi} are zero, because y is restricted to 0 < y < 1. Integrals over larger values of x simply multiply the answers just given by factors of Exp[-2 Pi n], n a positive integer. Hence, the complete answer is

Sum[Exp[-2 n Pi], {n, 0, Infinity}] (int1 + int2) //Simplify
(* (E^(2*Pi - ArcSin[y]) + E^(Pi + ArcSin[y]))/((-1 + E^(2*Pi))*Sqrt[1 - y^2]) *)

More Idiosyncrasies

In general, Integrate works correctly for any domain not including branch point of ArcSin[x]; e,g., for

Integrate[Exp[-x]*DiracDelta[y - Sin[x]], {x, (n - 1/2) Pi, (n + 1/2) Pi}, 
 Assumptions -> {0 < y < 1}]

with n any explicitly given integer. (However, if n is left an unspecified integer, n \[Element] Integers, Integrate returns unevaluated.) On the other hand, Integrate applied to domains including branch points returns unevaluated, with the exception of a domain extending to Infinity.

One might expect that replacing Sin[x] by Cos[x] would produce similar results. Indeed, Integrate applied to domains not including branch points of ArcCos[x] works fine, and applied to domains including branch points returns unevaluated. However, domains extending to Infinity also return unevaluated. So, the error found by @Markus occurs for Sin[x] but not for Cos[x].

@Jens observed in a Comment that Integrate seems to work well for domains including branch points, when the function in question is a polynomial. For instance, I have found that

Integrate[Exp[-x]*DiracDelta[y - (x - 1)^2], {x, 0, 2}, Assumptions -> {0 < y < 1}] 

returns the correct answer despite a branch point in the inverse of (x - 1)^2 at x = 1.

(* Cosh[Sqrt[y]]/(E*Sqrt[y]) *)

I also experimented with using BesselJ[0, x] and BesselJ[1, x], which oscillate somewhat like Cos[x] and Sin[x], respectively. However, Integrate returned unevaluated for all cases I tried.