Re: Integrating "If"
- To: mathgroup at smc.vnet.net
- Subject: [mg84531] Re: Integrating "If"
- From: "David W.Cantrell" <DWCantrell at sigmaxi.net>
- Date: Wed, 2 Jan 2008 06:57:36 -0500 (EST)
- References: <flc7vh$809$1@smc.vnet.net>
hemphill at alumni.caltech.edu wrote:
> I've come across what appears to be a bug in my (rather old) version of
> Mathematica, so if it's fixed in a newer version, please let me know.
I'm surprised that you haven't gotten a response before now.
> In[1]:= $Version
>
> Out[1]= 5.1 for Linux (October 25, 2004)
>
> This is a simplified version of the problem I am actually considering:
> Two points are chosen randomly from U(0,1). What is the distribution
> of d, the distance between them? If I formulate the problem this way,
> I get a correct answer:
>
> In[2]:= Integrate[If[Abs[x0 - x1] < d, 1, 0], {x0, 0, 1}, {x1, 0, 1},
> Assumptions -> {d > 0}]
>
> Out[2]= Piecewise[{{3/4, d == 1/2}, {1, d >= 1}}, 2*d - d^2]
BTW, it's rather unfortunate that the case d == 1/2 is considered
separately. It would have been nice if FullSimplify[Out[2]] were to give
Piecewise[{{1, d >= 1}}, 2*d - d^2] or Piecewise[{{1, d > 1}}, 2*d - d^2].
> I am delighted that Mathematica can integrate the "If" function, but
> if I use this equivalent (I think!) formulation, I get an incorrect
> answer:
>
> In[3]:= Integrate[If[(x0 - x1)^2 < d^2, 1, 0], {x0, 0, 1}, {x1, 0, 1},
> Assumptions -> {d > 0}]
>
> Out[3]=
> Piecewise[{{5/8, 2*d == 1}, {1, d >= 1}, {((4 - 3*d)*d)/2,
> 1/2 < d < 1}}, (-1 + 6*d - 3*d^2)/2]
>
> In[4]:= Plot[%, {d, 0, 2}, PlotRange -> All]
>
> This function of d should be positive, monotonic and continuous, but
> this plot exhibits none of those traits. I would prefer to use the
> second formulation, because it is easy to generalize for sampling from
> the unit square, unit cube, etc.
Version 6 gives the _same_ result for both your In[2] and your In[3],
namely, a correct result: your Out[2].
> By the way, I included the "Assumptions -> {d > 0}" in my formulations
> so that Mathematica would be aware that d was a real variable, and
> reduce the number of cases it had to consider.
Good idea.
> If I remove it from
> In[3], then I get correct results for positive d. However, the
> results for negative d, which should be a mirror image of those for
> positive d, are screwy.
In version 6, when Assumptions -> {d > 0} is not given, the result is
correct:
In[4]:= Integrate[If[(x0 - x1)^2 < d^2, 1, 0], {x0, 0, 1}, {x1, 0, 1}]
Out[4]= Piecewise[{{3/4, d == -(1/2) || d == 1/2},
{1, d >= 1 || d <= -1}, {-d^2 + 2*Sqrt[d^2],
-1 < d < -(1/2) || -(1/2) < d < 0 ||
0 < d < 1/2 || 1/2 < d < 1}}]
> Is there some way I can reformulate the integral, playing to
> Mathematica's strengths, so that I can get an answer I can trust?
I'm not sure. But you should be aware that there are alternative ways which
do not involve If. Both Boole and UnitStep exist in your version 5.1, and
so, equivalent to my In[4], you could try using either
Integrate[Boole[(x0 - x1)^2 < d^2], {x0, 0, 1}, {x1, 0, 1}]
or
Integrate[UnitStep[d^2 - (x0 - x1)^2], {x0, 0, 1}, {x1, 0, 1}]
In version 6, they both yield my Out[4].
> Keep in mind that I am also interested in evaluating integrals like:
>
> Integrate[If[(x0-x1)^2+(y0-y1)^2<d^2,1,0],{x0,0,1},{x1,0,1},{y0,0,1},{y1,
> 0,1}]
>
> (This integral takes a while to evalulate, and the result isn't
> continuous or monotonic. I'm about half way done with evaluating it
> by hand--it's an interesting problem.)
Yes, it does take a while, and I didn't have the patience to wait. But at
least you can see the result fairly quickly using, say,
Plot[NIntegrate[UnitStep[d^2 - (x0 - x1)^2 - (y0 - y1)^2], {x0, 0, 1},
{x1, 0, 1}, {y0, 0, 1}, {y1, 0, 1}], {d, 0, 2}]
BTW, what did you get by hand eventually?
David