1
$\begingroup$

I have an algorithm which needs to calculate the difference of two Gamma distributions evaluated at some large values. I cannot do this over CDFs because they give $1$ if a large value is evaluated. Then the difference of the CDFs is always $0$. Getting a difference of zero is breaking down my code at the optimzation stage. I am able to solve this problem if I use NIntegrate over the PDFs but every iteration needs alot of time and eventually the code is not running due to impracticality or in other words I must wait a day which normally should end in a few minutes at max. Here is my code:

m = 1;
ω = {10, 10, 10, 10, 10, 10, 10, 10, 10, 10};
γ = {0.1202, 0.1413, 0.1862, 0.2399, 0.3090, 1.0000, 1.5849, 1.9055, 1.9953, 2.5704};

f0[k_,y_] := PDF[ChiSquareDistribution[ω[[k]]], y]
f00[k_, z_] := Gamma[ω[[k]]/2, 0, z/2]/Gamma[ω[[k]]/2]

f1[k_,y_] := \[Piecewise]   ((2(γ[[k]]+1))^(-ω[[k]]/2) E^(-y/(2(γ[[k]]+1))) y^(-1+ω[[k]]/2))/Gamma[ω[[k]]/2]    y>0  0 True
f11[k_, z_] := Gamma[ω[[k]]/2, 0, z/(2 (γ[[k]] + 1))]/Gamma[ω[[k]]/2]

t={120,200};
p0[k_, idx_] := f00[k, t[[idx+1]]] - f00[k, t[[idx]]]
p1[k_, idx_] := f11[k, t[[idx + 1]]] - f11[k, t[[idx]]]
utmp[k_, idx_] := Log[p1[k, idx]/p0[k, idx]]
Subscript[m, 0][k_] := Sum[p0[k, idx]*utmp[k, idx], {idx, 1, m}];
Subscript[s, 0] := Sum[Sum[p0[k, idx]*(utmp[k, idx] - Subscript[m, 0][k])^2, {idx, 1, m}], {k, 1, 10}];

Lets evaluate the code:

Subscript[s, 0]

Indeterminate

If we change

t = {20, 25};

we get

0.517914

If I change

p0[k_, idx_] := NIntegrate[f0[k, y], {y, t[[idx]], t[[idx + 1]]}]
p1[k_, idx_] := NIntegrate[f1[k, y], {y, t[[idx]], t[[idx + 1]]}]

Then for the original case

t={120,200}; 

we get

3.17318*10^-17

But in this case as I said before evaluating Subscript[s, 0] is taking too much time.. Is there anyway to do with without NIntegrate?

here it isenter image description here

$\endgroup$
6
  • $\begingroup$ I'm having a hard time understanding the issue. The title is clear about the CDF of a Gamma distribution. But nowhere in the body is the use of the CDF function on a gamma distribution. The CDF of the Gamma distribution is exact if you use exact numbers: CDF[GammaDistribution[1, 3], 500] - CDF[GammaDistribution[1, 3], 499] gets one $\frac{1}{e^{499/3}}-\frac{1}{e^{500/3}}$. Using N on that expression gets one 1.64005*10^-73. $\endgroup$
    – JimB
    Feb 13 '20 at 16:30
  • $\begingroup$ @JimB f00 is the CDF of the density f0 and f11 is the CDF of the density f1. I am talking about the inaccuary of calculating p0[k_, idx_] if I use directly the CDF, namely the definition of f00 given in the code above. That is actually CDF of chi squared density not Gamma density as you mentioned but I just wrote it in the title having Gamma everywhere so the CDF of Chi-squared. f1 and f11 are also some kind of chi squared but I dont know exact details. $\endgroup$ Feb 13 '20 at 16:42
  • $\begingroup$ f0 is f0 and not CDF. Please show why using CDF[whatever distribution] doesn't work. Are you saying that your use of CDF in a function you created doesn't work or that the CDF function doesn't work as explicitly stated in your title? $\endgroup$
    – JimB
    Feb 13 '20 at 16:51
  • $\begingroup$ @JimB I didnt say f0 is a CDF, f0 is a density function. f00 is a CDF. If you run my code you will see that p0[k_, idx_] := f00[k, t[[idx+1]]] - f00[k, t[[idx]]] is evaluated to zero. This result is totally wrong. The true result can be obtained by using p0[k_, idx_] := NIntegrate[f0[k, y], {y, t[[idx]], t[[idx + 1]]}] which is another way of calculating the same thing. I guess you are unhappy with the mismatch between the title and the question? $\endgroup$ Feb 13 '20 at 16:55
  • $\begingroup$ Not unhappy, just confused. There is a lot of code to wade through with no comments inserted AND the title doesn't match the code. Your code for f0 might be mathematically equivalent to CDF but CDF has a whole lot more going for it under the hood. Have you tried using the WorkingPrecision option? $\endgroup$
    – JimB
    Feb 13 '20 at 17:02
3
$\begingroup$
m = 1;
ω = {10, 10, 10, 10, 10, 10, 10, 10, 10, 10};

Convert γ to exact values to avoid forcing machine precision calculations

γ = {0.1202, 0.1413, 0.1862, 0.2399, 0.3090, 1.0000, 1.5849, 1.9055, 
    1.9953, 2.5704} // Rationalize;

f0[k_, y_] := PDF[ChiSquareDistribution[ω[[k]]], y]
f00[k_, z_] := Gamma[ω[[k]]/2, 0, z/2]/Gamma[ω[[k]]/2]

f1[k_, y_] := 
   Piecewise[{{y^(-1 + ω[[k]]/2)/
             ((2*(γ[[k]] + 1))^(ω[[k]]/2)*
                E^(y/(2*(γ[[k]] + 1))))/
          Gamma[ω[[k]]/2], y > 0}}]

f11[k_, z_] := 
 Gamma[ω[[k]]/2, 0, z/(2 (γ[[k]] + 1))]/Gamma[ω[[k]]/2]

t = {120, 200};
p0[k_, idx_] := f00[k, t[[idx + 1]]] - f00[k, t[[idx]]]
p1[k_, idx_] := f11[k, t[[idx + 1]]] - f11[k, t[[idx]]]
utmp[k_, idx_] := Log[p1[k, idx]/p0[k, idx]]
Subscript[m, 0][k_] := Sum[p0[k, idx]*utmp[k, idx], {idx, 1, m}];
Subscript[s, 0] := 
  Sum[Sum[p0[k, idx]*(utmp[k, idx] - Subscript[m, 0][k])^2, {idx, 1, m}], {k, 
    1, 10}];

In evaluating Subscript[s,0] use arbitrary-precision rather than machine precision

Subscript[s, 0] // N[#, 15] &

(* 3.17318496695314*10^-17 *)

EDIT : Re your comment

If you enter machine precision numbers you will get machine precision results.

Precision[104.93434075028202]

(* MachinePrecision *)

N[f00[9, 200] - f00[9, 104.93434075028202], 100]

(* 0. *)

Precision[%]

(* MachinePrecision *)

Use either

f00[9, 200] - f00[9, 104.93434075028202`40]

(* 5.5830137827358486462429*10^-18 *)

Precision[%]

(* 22.7469 *)

Note that the complexity of the calculation results in a loss of about 17 digits of precision

or

f00[9, 200] - f00[9, SetPrecision[104.93434075028202, 40]]

(* 5.5830137827358526192441*10^-18 *)

Precision[%]

(* 22.7469 *)

EDIT 2: I cannot reproduce the results that you show in your comments. I am using version 12.0.0 on a Mac.

f11[1, SetPrecision[525/4, 40]] - f11[1, SetPrecision[225/2, 40]]

(* 4.4711124996345361867775*10^-17 *)

Precision[%]

(* 23.3494 *)

N[f11[1, 525/4] - f11[1, 225/2], 30]

(* 4.47111249963453618677749651064*10^-17 *)

Precision[%]

(* 30. *)
$\endgroup$
12
  • $\begingroup$ It works but not always. Thats really very very strange.. I received an error message because this N[f00[9, 200] - f00[9, 104.93434075028202], 100] was evaluated to zero although NIntegrate[f0[9, y], {y, 104.93434075028202, 200}] gives me 5.58301*10^-18 $\endgroup$ Feb 14 '20 at 10:50
  • $\begingroup$ do you have any idea what is special with this number 104.93434075028202? if I replace this number with 105 again I get something which is small but non-zero so works fine but not for this special number.. $\endgroup$ Feb 14 '20 at 10:52
  • $\begingroup$ Bob, thank you very much for the edit. You are right but It seems that it is out of my control because of my lack of mathematica knowledge. I am using NMaximize and inside this, there is f00[...]-f00[x] and here, $x$ is the variable over which NMaximize is maximizing. How can I tell NMaximize to stop using machine precision but rather to use the precision of my interest? else it finds at some point of maximization some indeterminate results and quits the optimization.. $\endgroup$ Feb 15 '20 at 0:27
  • $\begingroup$ I did like this NMaximize[{r[k, idx, [Lambda]], 0 <= [Lambda] <= 149.9}, [Lambda], WorkingPrecision -> 30] and i still have the same problem but it says "The precision of the argument function (0<=[Lambda]<=149.9) is less \ than WorkingPrecision (30)." and NMaximize::nnum: The function value Indeterminate is not a number at {[Lambda]} = {0.170947490850151928576394766424}. $\endgroup$ Feb 15 '20 at 17:29
  • $\begingroup$ 030 <= \[Lambda] <= 149.930 after this I dont get this one "The precision of the argument function (0<=[Lambda]<=149.9) is less \ than WorkingPrecision (30)." but I still get The function value Indeterminate is not a number at {[Lambda]} = {0.170947490850151928576394766424} $\endgroup$ Feb 15 '20 at 17:35
1
$\begingroup$

Edit: I've changed this from an extended comment to an answer.

The title concerns estimating the difference between two $\chi^2$ distribution functions. The text has the appearance of not directly addressing that question because of the large amount of code that doesn't explicitly mention a distribution function (at least in my opinion).

A brute force approach is to create a function that calculates the difference explicitly:

f[n_?IntegerQ, z1_?NumericQ, z2_?NumericQ, nDigits_: 50] := 
 N[CDF[ChiSquareDistribution[n], Rationalize[z2]] - 
   CDF[ChiSquareDistribution[n], Rationalize[z1]], nDigits]

For an example:

f[10, 199, 200]
(* 9.9470513771808598331974744855157686475664682149287*10^-38 *)

However, this direct approach fails for more extreme values of $z_1$ and $z_2$. Note that the difference in the CDF functions is a function of the incomplete gamma function:

PiecewiseExpand[
   CDF[ChiSquareDistribution[n], z2] - CDF[ChiSquareDistribution[n], z1],
   Assumptions -> {z1 > 0, z2 > 0}] // FunctionExpand // FullSimplify


(* (Gamma[n/2, z1/2] - Gamma[n/2, z2/2])/Gamma[n/2] *)

This is the same as using

Gamma[n/2, z1/2, z2/2]/Gamma[n/2]

So a more stable function for the difference in two $\chi^2$ distribution functions (with the same number of degrees of freedom) is

g[n_?IntegerQ, z1_?NumericQ, z2_?NumericQ, nDigits_: 50] := 
 N[Gamma[n/2, Rationalize[z1/2], Rationalize[z2/2]]/Gamma[n/2], nDigits]
g[10, 199999, 200000]

(* 9.6305672408760073406090960801307527696848945201455*10^-43412 *)
$\endgroup$
4
  • $\begingroup$ Hello Jim. Would you mind trying the code f11[1, SetPrecision[525/4, 40]] - f11[1, SetPrecision[225/2, 40]] at your mathematica and let me know if it gives zero or something like 10^-17? The issue is about the answer posted by Bob. $\endgroup$ Feb 16 '20 at 0:50
  • $\begingroup$ I get 4.4711124996345361867775*10^-17 on Windows 10, Mathematica 12.0. $\endgroup$
    – JimB
    Feb 16 '20 at 3:35
  • $\begingroup$ It would still be good to make an explicit match with your title and your code. Those that might be most helpful to you don't need to know the form of a chisquare cumulative distribution function. In other words, make an minimal working example. I still don't see that you've done that. You essentially have "CDF" in your title but the code only has "PDF". $\endgroup$
    – JimB
    Feb 16 '20 at 3:45
  • $\begingroup$ Jim, thank you very much for the answer. I also got the same result with copy paste but with my own code I still get 0 and both are actually the same.. anyways.. ?Numeric is a very good idea.. $\endgroup$ Feb 16 '20 at 20:34

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.