Re: Integrating "If"
- To: mathgroup at smc.vnet.net
- Subject: [mg84543] Re: Integrating "If"
- From: Scott Hemphill <hemphill at hemphills.net>
- Date: Thu, 3 Jan 2008 05:35:52 -0500 (EST)
- References: <flc7vh$809$1@smc.vnet.net> <flfubb$jrm$1@smc.vnet.net>
- Reply-to: hemphill at alumni.caltech.edu
"David W.Cantrell" <DWCantrell at sigmaxi.net> writes:
> 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].
For me, both of these methods yield the same buggy result as using "If".
>> 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}]
This didn't work for me. The integrator complained that the result
coverged too slowly.
> BTW, what did you get by hand eventually?
If 0 <= d <= 1, then
(-8*d^3)/3 + d^4/2 + d^2*Pi
If 1 <= d <= Sqrt[2], then
1/3 - 2*d^2 - d^4/2 + (4*Sqrt[-1 + d^2]*(1 + 2*d^2))/3 + Pi +
2*(-1 + d^2)*ArcCot[Sqrt[-1 + d^2]] - 2*(1 + d^2)*ArcTan[Sqrt[-1 + d^2]]
I still used Mathematica to do some of the algebra. I had to coerce
Mathematica to give the answer in this form, since it wanted to return
answers with "I" and "ArcSinh" in them.
I had started to use a complicated method which evaluated areas of
circles which were clipped by the unit square. Maxim Rytin (private
corresponence) pointed out that the pdf of (x0-x1)^2 is (1/Sqrt[u]-1),
0<u<=1, and to get the pdf of (x0-x1)^2+(y0-y1)^2, you merely have to
convolve the first pdf with itself. After that, you integrate to get
the cdf, and substitute d^2 for u. A little messy, but not too bad.
Scott
--
Scott Hemphill hemphill at alumni.caltech.edu
"This isn't flying. This is falling, with style." -- Buzz Lightyear