hanchenphd/GMMAT

Saddlepoint approximation in `.pKuonen()`

Closed this issue · 4 comments

dwoll commented

First of all, many thanks for making your work available for others to use and study! I'm trying to study implementations of the saddlepoint approximation from Kuonen (1999), and I'd be glad if you could help me with the following questions.

In function .pKuonen(), it seems to me that the following line does not achieve anything. Maybe a return() is missing here?

GMMAT/R/SMMAT.R

Line 1170 in 518f664

pchisq(x/lambda, df = df, ncp = delta, lower.tail = FALSE)

In sub-functions kprime0()

GMMAT/R/SMMAT.R

Line 1181 in 518f664

sum(((delta + df) * lambda)/(1 - 2 * zz * lambda) + 2 * (delta *

and kpprime0()

GMMAT/R/SMMAT.R

Line 1186 in 518f664

sum((2 * (2 * delta + df) * lambda^2)/(1 - 2 * zeta * lambda)^2 + 8 *

The code is using a different notation compared to Kuonen (1999) p. 930. Kuonen uses sigma^2 in the equations (non-centrality parameter for a chi^2 distribution) whereas you use df+delta. My math skills are quite limited, so it would help me a lot if you could explain the relationship between your implementation and Kuonen's equations. Many thanks in advance!

Thanks for catching line 1170. Yes, this line of code should be wrapped by return() and we will fix it in the next version.

The kprime0 and kpprime0 functions are the first and second derivatives of the cumulant generating function K(ζ). Kuonen used h_i and \sigma_i^2 to denote the degree of freedom and the non-centrality parameter, and they were df and delta in the code. I think there was a missing ζ in the numerator of the second term in Kuonen's K(ζ), and perhaps that could explain the confusion when you take the derivatives.

Best,
Han

dwoll commented

Thank you very much for your quick response and for taking the time to help me better understand your implementation!

In Kuonen's PhD thesis equation 3.2.4 (page 21), ζ is also missing from the numerator in the second term of the cumulant generating function K. The first and second derivatives on page 24 also look somewhat different compared to your implementation. That is indeed what is confusing to me. Is there a manuscript where the formulas you are implementing are written down? Thanks again!

I believe it was a typo in Kuonen's paper. The cumulant generating function is basically the log of the MGF, and the MGF of a weighted sum of independent non-central chi-squares is the product of the MGF of each term. See Wikipedia and StackExchange for the MGF of non-central chi-squares.

Kuonen cited Equation 2.3 from Imhof 1961, which is consistent with the links above.

Best,
Han

dwoll commented

That makes sense, thank you for the explanation!