3

Motivation/background: The following integrals were part of an argument in my PhD thesis that a certain piece of forensic evidence has no inculpatory value. I gave up trying to solve them analytically using Mathematica, and settled instead for numerical estimates without error bounds. Now I want to turn the argument into a computer-checkable proof, from a mix of objective and subjective assumptions. Since it is a complicated argument, I would like to be able to minimize the number of assumptions as much as possible. For that, I need either an analytic solution, or a numerical solution that includes provable error bounds. If there is a better place for this, I would be grateful to know it (I have other problems that I cannot solve, and will generate more in the future).

Actual question: The integrals are of the following form, for a few different settings of the parameter $S$. I expect that anyone who can solve it for one such $S$ can solve it for the others, so I'll give just one example: $S$ is the set of $$\langle p_1,p_2,p_3 \rangle \text{ such that } 0 < 1-p_1-p_2-p_3 < p_3 < p_2 < 1$$ Then I want to know how to solve analytically, or numerically with a provable error bound, this: $$ \int_{\alpha \in (.85,1)} \int_{\langle p_1,p_2,p_3\rangle \in S} f(p_1,p_2,p_3,\alpha) \ dp_1 \, dp_2 \, dp_3 \, d\alpha $$ where $f(p_1,p_2,p_3,\alpha)$ is: $$ (p_1\alpha + \frac{89}{90}(1-\alpha))^{10}\ (p_2\alpha + \frac{1}{90}(1-\alpha))^{20}\ (p_3\alpha)^{40} \ ((1-p_1-p_2-p_3)\alpha)^{19} $$

I would also appreciate suggestions for making it easier to solve by a computer algebra system running on a university server or cluster.

Dustin Wehr
  • 285
  • Try Math.Stackexchange. This is not the right place for such a problem. – Ben McKay Apr 10 '16 at 19:25
  • @BenMcKay Perhaps I'm just betraying my ignorance here, but I'd be interested in seeing suggestions/solutions from people in numerical-analysis-influenced fields, since I've got some similarly annoying integrals I'd like to estimate as background to some of my own research. (The high exponents made my naive attempts to crunch the integrals with a CAS prohibitively slow) –  Apr 10 '16 at 19:31
  • 1
    This is just a polynomial, right? An analytic solution should be trivial. – Nate Eldredge Apr 10 '16 at 19:31
  • 3
    Maybe I am being to hasty. You only have a 10th power of something messy if you work with the variable $q_1=1-p_1-p_2-p_3$ instead of $p_1$, instead of a 19th power, and the domain is a little easier. – Ben McKay Apr 10 '16 at 19:34
  • 3
    @NateEldredge Wouldn't the same logic say that integrating polynomial expressions of matrix coefficients of e.g. SU(2) should be trivial? (And if that is trivial, then my next question is: can I send you an email with some integrals I'd like estimated... ) –  Apr 10 '16 at 19:35
  • @YemonChoi: In that case it's made difficult by the measure. Here the measure is Lebesgue. And by "trivial" I of course mean "for a computer". I don't understand why the OP's Mathematica attempts failed. – Nate Eldredge Apr 10 '16 at 19:36
  • 2
    @NateEldredge Well by identifying SU(2) with the 3-sphere in R^4 and using homogeneity (+ a trick with a Gaussian weight) I can convert to an integral over R^4, but it doesn't make the polynomials involved much nicer, especially since I need lower bounds rather than upper bounds. (But this is getting off-topic -- my contention was merely that "calculable in principle using 1st-year calculus" doesn't always translate into "calculable in practice in an efficient re-usable way".) –  Apr 10 '16 at 19:39
  • @YemonChoi In this particular example, Maple seems to be both efficient and reusable. Of course, we are not told what $S$ is; it may screw things up. –  Apr 10 '16 at 20:01
  • Thank you all. I will try on math.stackexchange. I also added a clarification to the end: If I could get a computer algebra system to return decent bounds on the answer, which are correct provided the system's code is correct, after hours or days of crunching, that would be enough.

    @AlexDegtyarev S is given in the first display-math line.

     – Dustin Wehr Apr 10 '16 at 20:02
  • 4
    In that particular example of $S$ I think that Maple will give the result within fractions of a second :) And there's no need to "code" it: you just write int() four times :) –  Apr 10 '16 at 20:08
  • @DustinWehr: Please don't post to Math.SE right away. This question can be migrated there. If you post a new question there, they will be duplicated, and people who see one will miss the progress made on the other. – Nate Eldredge Apr 10 '16 at 20:13
  • I would love to hear how estimating this integral implies that "a certain piece of forensic evidence has no inculpatory value"... –  Apr 10 '16 at 20:21
  • sage gives 0.00014320685701723506 in seconds. I can put up the sage code if you like, but I didn't do anything clever. – Ben McKay Apr 10 '16 at 20:25
  • 2
    @MarkLewko it's pages 58-73 of http://arxiv.org/pdf/1502.02272v1.pdf – Dustin Wehr Apr 10 '16 at 20:33
  • 2
    @YemonChoi Perhaps you should post your $SU(2)$ integrals! Ordinarily, I would think that an individual polynomial of the type you describe shouldn't be bad to integrate exactly, though asymptotics of a family might be hard to obtain. – David E Speyer Apr 11 '16 at 00:05

2 Answers2

6

Using the change of variables suggested by Ben McKay, I did the following in Mathematica:

f[q1_, p2_, p3_, a_] = (-(q1 - 1 + p2 + p3)*a + (89/90)*(1 - a))^10*(p2*
     a + (1/90)*(1 - a))^20*(p3*a)^40 *(q1*a)^19

f2[q1_, p2_, p3_] = Integrate[f[q1, p2, p3, a], {a, 85/100, 1}]
f3[p2_, p3_] = Integrate[f2[q1, p2, p3], {q1, 0, p3}]
f4[p2_] = Integrate[f3[p2, p3], {p3, 0, p2}]
ans = Integrate[f4[p2], {p2, 0, 1}]

The result is exactly

33672470631114496923827791311023515239260091504681177509643881461349318033718604944405779863067208073100207105031683661728771287430823858857293983479827827539912836287902935668145188731720876839417022184157831269987901 / 537416971039477383994538150480683103634950206009134466653124329302866504399642121590258750479296102400000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

which is, to within 50 decimal places,

0.000062656135637073129567523815922285464561331045177648

This took less than a minute on an ordinary laptop.

Nate Eldredge
  • 97,710
3

Edit: (11 April 7:03 est)

In sage

a=var('a')
q = list(var('q_%d' % i) for i in range(1,4))
f=((1-(q_1+q_2+q_3))*a+(89/90)*(1-a))^(10)*(q_3*a+(1/90)*(1-a))^(20)*(q_2*a)^40*(q_1*a)^(19)
integral(integral(integral(integral(f,q_1,0,q_2),q_2,0,q_3),q_3,0,1),a,.85,1)

yields $6.265613563707314 \times 10^{-05}$

Ben McKay
  • 1,128
  • Great! You said you did nothing clever, but was the change of variables necessary in Sage? Would you mind briefly explaining why you thought it would make such a difference? – Dustin Wehr Apr 10 '16 at 20:44
  • 1
    As pointed out by Robert Israel, the region should be $0 < q_1 < q_3 < q_2 < 1$ (unless that was a typo by the OP). Also, your definition of $f$ has $q_3$ twice; the first one should be $q_2$. Making those changes, I get the same result in Sage as in my Mathematica below. See https://cloud.sagemath.com/projects/17eb8257-e5ad-41b8-889a-fce73e2a41e9/files/mo%20integral.sagews – Nate Eldredge Apr 11 '16 at 02:53
  • @NateEldredge: right. I fixed it above, and now it agrees with your ansswer. – Ben McKay Apr 11 '16 at 06:06
  • 2
    @DustinWehr: The change of variables might not have been necessary. The point of it is to look at the highest terms and try to make them pure powers of single variables, so that when you expand out (using the binomial theorem) each of the terms in $f$ you don't get so many terms in the resulting polynomial. – Ben McKay Apr 11 '16 at 12:49