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