12
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. $\endgroup$ – bbgodfrey Jan 18 '15 at 5:43
  • 4
    $\begingroup$ This seems to be a bug. As you correctly point out, there are infinitely many zeros in the argument of the delta function, and they must be accounted for in a sum. There is no justification for omitting that sum, and it would give you an infinite series. $\endgroup$ – Jens Jan 18 '15 at 6:45
  • 1
    $\begingroup$ I've reported this as a bug (CASE:2323185). $\endgroup$ – Jens Jan 19 '15 at 0:26
  • $\begingroup$ The result in 10.4.1 is still wrong, but at least different. $\endgroup$ – J. M. will be back soon Jul 10 '16 at 18:21
  • $\begingroup$ @J.M. Strange - in version 11 (prerelease), I don't get anything different - still the same wrong result... $\endgroup$ – Jens Jul 10 '16 at 22:18
11
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ (+1) Interestingly, the result looks correct if one does it with a polynomial instead of Sin[x], so the existence of branches of the inverse is not per se a deal-breaker. $\endgroup$ – Jens Jan 18 '15 at 17:40
  • $\begingroup$ Yes, @Jens. I see that (x-1)^2 instead of Sin[x] appears to work. Perhaps, Integrate can handle branch points in polynomials but not in transcendental function. In any case, someone should report this to Wolfram, Inc. $\endgroup$ – bbgodfrey Jan 18 '15 at 23:43
  • $\begingroup$ Disturbingly, the functions in the integral do not have a branch cut. Inverting the sine causes the problem, and so one actually needs to know how to evaluate such integrals in order to spot the trouble. Why doesn’t Mathematica at least throw a warning flag at me, as it occasionally does when branch cut issues arise? I’d be happy to forward this to Wolfram Inc., but have never done such a thing. Do they have some standardized bug-report procedure? $\endgroup$ – Markus Deserno Jan 19 '15 at 0:29
  • $\begingroup$ @Markus, I suggest you contact Wolfram on its web page. Do not be surprised, it you do not receive a response for a few days. Include a reference your StackExchange question. $\endgroup$ – bbgodfrey Jan 19 '15 at 17:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.