3
$\begingroup$

Similar things have been asked previously (see https://academia.stackexchange.com/q/71536/4132 for instance), but I want to address the specific case in which a complex symbolic/analytical result has been returned by Mathematica. How can I check, verify or justify that such a result is right, in the context of scientific research?

This is my concrete case: I was interested in calculating the expected value of $X^{-1}$ and $X^{-2}$, where $X$ follows a noncentral $\chi^2$ distribution. I tried to get an analytical or symbolic result with Mathematica. This is what I did:

Xm1[nu_, delta_] := 
 TransformedDistribution[x^(-1), x \[Distributed] X[nu, delta]]

Mean[Xm1[nu, delta]]

After a while, I fortunately got this result:

-2^(-2 + nu/2) delta^(
 1 - nu/2) E^(-(delta/2) + (I nu \[Pi])/
  2) (Gamma[-1 + nu/2] - Gamma[-1 + nu/2, -(delta/2)])

(which involves the use of complex numbers!)

Something similar was done to obtain an expression for $E[X^{-2}]$.

I've numerically checked (for concrete values of $\nu$ and $\delta$) that the result returned by Mathematica seems to be correct (or, at least, consistent with other computations done with Mathematica).

I would like to be able to justify this analytical result in case I use it in a scientific paper. I am not sure whether it is enough to say: "This was obtained by intense symbolic computations [done with Mathematica]".

What else can I say to prove that the analytical expression for $E[X^{-1}]$ is right?


UPDATE

My question is quite generic, but I am actually interested in the case I've shown in the example above,that is, $X \sim \chi^2_{\nu}(\delta)$ ($X$ follows a noncentral $\chi^2$ with $\nu$ degrees of freedom and noncentrality parameter $\delta$).

$\endgroup$
  • $\begingroup$ Can you give a complete example please (i.e. define X)? $\endgroup$ – Szabolcs Sep 13 '17 at 10:21
  • $\begingroup$ @Szabolcs That's what I tried to do when I speak about a noncentral $\chi^2$ distribution. $\endgroup$ – Vicent Sep 13 '17 at 11:30
6
$\begingroup$

You're asking how to verify a definite integration? As far as I know, the only way to verify definite integration (you've thrown away information) is to redo it in a way that you trust. If you don't trust Mathematica's symbolic integration, you can always try to do it by hand, use a different solver, or integrate various controlled approximations to it.

Let's see what you're doing. Translating your text into an expression:

Xm1[nu_, delta_] := 
 TransformedDistribution[x^(-1), 
  x \[Distributed] NoncentralChiSquareDistribution[nu, delta]]

To calculate the Mean we're doing something like:

Integrate[ x PDF[Xm1[nu, delta]][x], {x, 0, \[Infinity]}, 
 Assumptions -> delta >= 0 && nu > 0]

Which indeed lands on your answer conditional upon nu>2.

Of course one way of disproving it is to have it fail consistency checks, like take $\nu$ and $\delta$ to regions where the mean should vanish or become singular. If your expression doesn't vanish / become singular in these limits, you should be suspicious.

As a verification game you can play using mathematica (or by hand): you can integrate a series expansion around vanishing delta to say third order:

approxIntegrand=Assuming[x > 0, 
 x * Normal[Series[PDF[Xm1[nu, delta]][x], {delta, 0, 3}]] // 
   ExpandAll // Collect[#, delta, FullSimplify] &]

(2^(-2 - nu/2) delta E^(-(1/2)/x) x^(-1 - nu/2) (1 - nu x))/ Gamma[1 + nu/2] + ( 2^(-5 - nu/2) delta^2 E^(-(1/2)/x) x^(-2 - nu/2) (1 + (2 + nu) x (-2 + nu x)))/Gamma[2 + nu/2] + ( 2^(-7 - nu/2) delta^3 E^(-(1/2)/x) x^(-3 - nu/2) (1 + (4 + nu) x (-3 + (2 + nu) x (3 - nu x))))/( 3 Gamma[3 + nu/2]) + (2^(-nu/2) E^(-(1/2)/x) x^(-nu/2))/Gamma[nu/2]

and match it against the same series expansion of your result.

Of course this isn't a proof -- it merely demonstrates consistency in a region you have control over.


Also I'm not sure if typing: Mean[Xm1[mu,delta]] precisely counts as intense symbolic computations ;-). But, as an aside, in Version 11.1 OS X at least, when I evaluate:

Assuming[delta >= 0 && nu > 2, Mean[Xm1[nu, delta]]]

I don't get your explicit complex dependence, instead it resolves happily to:

-2^(-2 + nu/2) (-delta)^(-nu/2) delta E^(-delta/2) (Gamma[-1 + nu/2] - Gamma[-1 + nu/2, -(delta/2)])


Update: to further discuss "how to justify" as per comment. I'm not so sure this adjunct is really a Mathematica answer as much as it is "good practices."

As a scientist you are justified to disseminate a result you believe in as much as you track and communicate any qualifications. Numerical support is sufficient for certain problems in certain fields. Professional symbolic software integration can generally be relied upon, but of course even Mathematica has had fairly spectacular bugs, so it's a good idea to cross-check and rederive in a number of ways to the point you're convinced.

My field is full of types that would be easily embarrassed when there's a nice closed form for an integration (in one variable!!) we couldn't eventually work out on our own (especially once we know the answer), even if the original integrand was intimidating and full of special functions. The nice thing is that almost everything intimidating is itself a solution to a differential equation, so if we know there's a closed path there's a hair to pull on. So we'll put extra energy into tracking it down -- relearning high-school math sometimes in the process. Then we'll put extra energy into making sure the end form is as compact/symmetric/etc as possible.

Often times however (and I'm thinking multi-loop Feynman integrals) we do rely on fairly sophisticated technology (like that of the Smirnovs ) and then we have to come up with our own consistency checks which we absolutely discuss and disclose.

The most important thing to do is convince yourself the expression is valid and simplified to a point you won't find embarrassing. Once you're there, depending upon pride/shame/dignity etc, it's possible to cite in following form:

"It is easy to see\footnote{e.g. using Mathematica, \tt{ Mean[xxx]}} that blah blah blah."

In general if you would feel comfortable trusting a table of definite integrals (do people feel comfortable trusting a table of definite integrals??) then you can cite as you would that.

If the only way you know how to get at a result is via code, you should make the algorithm clear, and for the Good of Science$^{\rm TM}$ consider sharing the code directly, and if you have any doubts, consider disclosing what you checked and how you convinced yourself that things were ok in the end. If the algorithm is a black box like Mathematica's symbolic integration then reference the black box. In general for 1D integrals, as a professional, it's probably worthwhile figuring out how to at least sketch out the integration -- convincing yourself it's true. When I get stuck, I ask my friends.

Aside: I had a collaborator with the following strategy: when he was blocking on something possibly silly (math/conceptual/etc) he would start by asking undergrads, then grad students, then postdocs, and eventually if the problem was tricky enough he was having brilliant conversations with his senior colleagues. This is not (merely) as cynically careerist as it might appear; it created a pretty thermal environment. Oftentimes people up and down would learn something they otherwise wouldn't have had he just sat in his room banging his head against the wall.

$\endgroup$
  • $\begingroup$ Thank you very much, @John! Some things to highlight: (1) My question is: How can I justify (in a scientific paper) that the expression obtained via Mathematica is true? I mean, I could try to deduce it by hand, but I do not know if this would be too hard. Is there a way to know the 'steps' that Mathematica has followed to reach to this result? I guess it won't be easy, since it takes several seconds to get the result (at least in my Windows 10 machine). (2) Thank you for the hint about getting rid of the complex part. $\endgroup$ – Vicent Sep 13 '17 at 13:04
  • $\begingroup$ And (3): About the expansion suggestion, thank you too, and: How did you get to manage Mathematica in such an efficient way?? $\endgroup$ – Vicent Sep 13 '17 at 13:04
  • $\begingroup$ @Vicent, "Is there a way to know the 'steps' that Mathematica has followed to reach to this result?" - no, unfortunately not, simply because software will not integrate functions in the same way a human does. If you're truly paranoid, don't trust possibly buggy software to do the algebra for you; do it by hand. $\endgroup$ – J. M. is away Sep 13 '17 at 13:25
  • $\begingroup$ @J.M. Not that I am paranoid, but I rather would like to give a mathematical reason for that result, other than "Mathematica says it". $\endgroup$ – Vicent Sep 13 '17 at 13:53
  • $\begingroup$ "a mathematical reason" - all the more reason to go through the algebra yourself, then. $\endgroup$ – J. M. is away Sep 13 '17 at 14:00
2
$\begingroup$

Making use of Maple 2017.2, I obtain $$E[X^{-1}]=1/2\,{\frac {{{\rm e}^{-\delta/2}} \left( \delta\, \left( -\delta \right) ^{-\nu/2}{2}^{\nu/2}\Gamma \left( \nu/2,-\delta/2 \right) - \delta\,\Gamma \left( \nu/2 \right) \left( -\delta \right) ^{-\nu/2}{ 2}^{\nu/2}+2\,{{\rm e}^{\delta/2}} \right) }{\nu-2}} $$ and $$E[X^{-2}]=1/4\,{\frac {{{\rm e}^{-\delta/2}} \left( {\delta}^{2}\Gamma \left( \nu/2 \right) \left( -\delta \right) ^{-\nu/2}{2}^{\nu/2}+\delta\, \Gamma \left( \nu/2 \right) \left( -\delta \right) ^{-\nu/2}{2}^{\nu/ 2}\nu-{\delta}^{2} \left( -\delta \right) ^{-\nu/2}{2}^{\nu/2}\Gamma \left( \nu/2,-\delta/2 \right) -\delta\, \left( -\delta \right) ^{- \nu/2}{2}^{\nu/2}\Gamma \left( \nu/2,-\delta/2 \right) \nu-4\,\delta\, \Gamma \left( \nu/2 \right) \left( -\delta \right) ^{-\nu/2}{2}^{\nu/ 2}+4\,\delta\, \left( -\delta \right) ^{-\nu/2}{2}^{\nu/2}\Gamma \left( \nu/2,-\delta/2 \right) -2\,{{\rm e}^{\delta/2}}\delta+4\,{ {\rm e}^{\delta/2}} \right) }{ \left( -4+\nu \right) \left( \nu-2 \right) }} $$

$\endgroup$
  • $\begingroup$ The result given by Maple is a little bit more 'intuitive'. Thank you. $\endgroup$ – Vicent Sep 13 '17 at 14:05
  • $\begingroup$ Does Maple give the result directly in LaTeX? If not, would you mind to include the raw output of Maple here? $\endgroup$ – Vicent Sep 13 '17 at 14:06
  • $\begingroup$ @Vicent: Here dropbox.com/s/whute9sj90d1typ/means.pdf?dl=0 it is. $\endgroup$ – user64494 Sep 13 '17 at 15:27
  • $\begingroup$ Thank you for the information! $\endgroup$ – Vicent Sep 14 '17 at 13:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.