Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2005
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2005

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Integrate gives wrong answer

  • To: mathgroup at smc.vnet.net
  • Subject: [mg57105] Re: Integrate gives wrong answer
  • From: Maxim <ab_def at prontomail.com>
  • Date: Mon, 16 May 2005 01:29:47 -0400 (EDT)
  • References: <d615eu$nnb$1@smc.vnet.net> <d66tf7$nnj$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

On Sun, 15 May 2005 07:24:23 +0000 (UTC), Maxim <ab_def at prontomail.com>  
wrote:

> On Fri, 13 May 2005 03:03:58 +0000 (UTC), Brian Rogers <brifry at gmail.com>
> wrote:
>
>> Hi, I have a problem on which Integrate misses badly. The code consits
>> of:
>>
>> fmin[x_, p_]:=[p(1-(x-a+e)/(2e))^(p-1)]/2e
>>
>> which is the pdf of the minimum of p draws from a uniform distribution
>> on [a-e,a+e], and
>>
>> Integrate[Abs[Min[mi,mc]-Min[mj,mc]]*fmin[mc,k]*fmin[mi,ni-k]*fmin[mj,nj-k],
>> {mc,a-e,a+e},{mi,a-e,a+e},{mj,a-e,a+e},Assumptions->{k\el\_Integers,
>> ni\el\_Integers,nj\el\_Integers,1<=k<=Min[ni,nj],0<a-e<a+e<1}]
>>
>> This produces the output:
>>
>> \!\(k\ \((\(a +
>>       e\)\/ni - \(2\ e\)\/\(1 + ni\) + \(a + e\)\/\(k - ni - nj\) +
>> \(2\ \
>> e\)\/\(1 - k + ni + nj\))\)\)
>>
>> which is not symmetric betwenn ni and nj, as it should be.
>>
>> Any ideas?
>> Thanks!
>> Brian Rogers
>> brifry at gmail.com
>>
>
> This can be done using the symmetries of the integral. First, it's
> independent of a, and second, the integral over mi < mj is the same as  
> the
> integral over mj < mi with ni and nj interchanged; also we can assume  
> that
> ni <= nj. The condition 0 < a - e < a + e < 1 is equivalent to 0 < e <  
> 1/2
> && e < a < 1 - e and we can choose a == 0:
>
> In[1]:=
> fmin[x_, p_] := p*(1 - (x - a + e)/(2*e))^(p - 1)/(2*e)
>
> Integrate[Abs[Min[mi, mc] - Min[mj, mc]]*
>        fmin[mc, k]*fmin[mi, ni - k]*fmin[mj, nj - k],
>      {mc, a - e, a + e},
>      Assumptions -> {a == 0, 1 <= k <= ni <= nj, 0 < e < 1/2,
>        a - e < mi < mj < a + e}] //
>    Simplify;
>
> Integrate[%, {mi, a - e, a + e}, {mj, mi, a + e},
>    Assumptions -> {a == 0, 1 <= k <= ni <= nj, 0 < e < 1/2}];
>
> % + (% /. {ni -> nj, nj -> ni}) // FullSimplify
>
> Out[4]=
> 2*e*(1/(1 + ni) + 1/(1 + nj) - 2/(1 - k + ni + nj))
>
> which seems to agree with the results of the numerical integration.
>
> Generally it's a good idea to specify the version of Mathematica you're
> using; the above will work in version 5.1 only. In 5.0 you can carry out
> the first integration step with the help of PiecewiseIntegrate from
> http://library.wolfram.com/infocenter/MathSource/5117/ .
>
> However, Mathematica 5.1 still has problems with this type of integrals:
>
> In[5]:=
> Assuming[e > 0 && p > 0,
>    Integrate[x*(e - 1 - x)^p, {x, -1 - e, -1 + e}]]
>
> Out[5]=
> 0
>
> which is incorrect. Integrating on 0 < e < 1 and 1 < e separately gives
> the correct answer.
>
> In[6]:=
> Integrate[e^(k - ni - nj)*(e - mi)^(-1 - k + ni)*
>        (e - mj)^(-1 + nj)*(mi - mj) +
>      e^(k - ni - nj)*(e - mi)^(-1 - k + ni)*(e - mj)^(-1 - k + nj)*
>        ((e - mi)^(1 + k) - (e - mj)^k),
>    {mi, -e, e},
>    Assumptions -> {1 < k < ni < nj, 0 < e < 1/2, -e < mj < e}]
>
> Out[6]=
> ComplexInfinity
>
> In fact the integral converges. A trickier issue is that if we perform  
> the
> integration steps in a different order or for different values of a, the
> result of integration may contain expressions of the form Beta[z, -ni, ni
> - 1]. For integer ni>0 (which is what we need) this is correct only as
> Limit[Beta[z, -n, n - 1], n -> ni].
>
> Maxim Rytin
> m.r at inbox.ru
>

A better way is

Block[{a = 0},
   Assuming[{1 <= k <= ni <= nj, 0 < e < 1/2,
       a - e < mi < mj < a + e},
     Integrate[Abs[Min[mi, mc] - Min[mj, mc]]*
           fmin[mc, k]*fmin[mi, ni - k]*fmin[mj, nj - k],
         {mc, a - e, a + e}] //
       Simplify]
];

Then the second integration step takes significantly less time.

Maxim Rytin
m.r at inbox.ru


  • Prev by Date: Re: Re: NIntegrate-FindRoot acting up in version 5.1
  • Next by Date: Re: Crossing of 3D functions
  • Previous by thread: Re: Integrate gives wrong answer
  • Next by thread: Re: GramSchmidt problem