Message Boards Message Boards

0
|
6603 Views
|
9 Replies
|
0 Total Likes
View groups...
Share
Share this post:

How to get accurate results of InverseCDF function in Mathematica 9?

I need to calculate values of the InverseCDF function of the noncentral Student T distribution. This is used for constructing statistical tolerance intervals. The values most commonly used are tabulated, but occasionally one needs values not found in the tables, or needs to be able to update the value when the parameters change.

The function to compute is: InverseCDF[NoncentralStudentTDistribution[s1-1,d],cl]

Note: Because of a symmetry in this function, one could equally well calculate: InverseCDF[NoncentralStudentTDistribution[s1-1,-d],1-cl]. The result would be the same. But this trick has not helped with the difficulties described below.

Here s1 is related to your sample size, but in some applications it won't be an integer (but it is positive). The order of s1 could be in the tens, hundreds, thousands, etc.

d is the non-centrality parameter that is related to the sample size and to the probability you are trying to cover (p). Normally d= Sqrt[s1]*InverseCDF[NormalDistribution[0,1],p], where p can be something like 0.9, 0.95, 0.99, etc. (less than 1). For these values of p, the InverseCDF of the normal distribution is usually between one and two. The calculation of the InverseCDF for the normal distribution is not the problem.

In some applications, the sample size used for calculating d will be different than s1 (say s2), but for this example I used the same.

cl is the confidence level one is trying to cover, something like 0.9, 0.95, 0.99, etc. (less than 1).

I was trying this using p=cl=0.95 and using the same sample size for d, so that the result is just a function of the number I use for s1 (I used integers in this example). I checked whether the results were correct by plugging the answer back into CDF and checking whether one gets 0.95 (I also used alternative checks to confirm that the errors I found are in the calculation of the InverseCDF, rather than in CDF itself). I noticed the following:

1) The calculation is fast and accurate for small sample sizes (<100).

2) For some values (e.g., 150) I get warnings about convergence issues in NIntegrate, but the result is accurate. The warning may go away for larger values (200). The calculation starts taking long (~ 1minute, although I did not time it accurately).

3) For larger values (420) I get inaccurate results (plugging into CDF gives ~0.503 instead of 0.95) with no warnings about convergence.

I tried adjusting the options in NIntegrate but it did not help.

I tried coding a simple bisection procedure to invert the CDF function and the results are much better. The calculation is faster and the results are accurate even when s1 is in the thousands, but I still get warnings about convergence. The difficulty in coding the bisection is finding an interval that always brackets the solution, but that is fairly tight. I found such bounds by looking into the theory of the non-central student distribution, but there is probably room for improvement.

I tried to clarify what is causing the numerical difficulties. The CDF is the integral of the PDF. As the sample size (s1) grows, the peak of the PDF is not getting narrower (I think I found that the half-width approaches a finite value), but the peak is moving farther and farther away from zero, so, it is becoming narrower relative to the intervals being considered in the integration.

Clearly this is a bug, since it gives an incorrect result without printing a warning. It is also clear that Mathematica could do the calculation better (as when I code a bisection procedure).

Is there a way to influence the way Mathematica inverts the CDF function?

Is there a different way to do this in the statistical packages (I am not very familiar with them myself)?

Another thing I plan to check is whether version 10 gives better results.

Thanks,

O.L.

POSTED BY: Otto Linsuain
9 Replies

Hello Jim,

Thanks very much for providing the article. I wrote to Wolfram, and they replied back. They thanked me and said they will work on the issue (contrary to some comments I have read on the web that they ignore bug reports). Of course, they have not fixed it yet.

I will read the article in detail and will probably use it. I found another article (which I have in my work computer, but which I can send if you are interested). That article provided two good approximations which I could then use as a starting point in FindRoot to get fairly good results. We need these values at work and what is published in tables is not always enough or is too inconvenient.

I looked into the issue myself and it is fairly technical. Some features of this distributions are known exactly, for others only bounds are known, but those could be used as a first guess or bracketing step in a solution method. I am not surprised that someone devoted a fundamental article just to the handling of this problem. I am a bit surprised that Wolfram let it slip. They could have been rushing through that development. Hopefully they will flush it out; they got smart people there.

Thanks a lot,

OL.

POSTED BY: Otto Linsuain
Posted 9 years ago
Posted 9 years ago

Excellent! That is a concise example that shows the issue. (One can format Mathematica code using the "Code Sample" icon. Just hover over the editing icons just below the "Reply to this discussion" text to see approximately what each means.) Here is your code doing so:

a = InverseCDF[NormalDistribution[0, 1], 0.95];
f[n_, cl_] := InverseCDF[NoncentralStudentTDistribution[n - 1, Sqrt[n]*a], cl]
fInv[n_, k_] := CDF[NoncentralStudentTDistribution[n - 1, Sqrt[n]*a], k]
Column[(AbsoluteTiming[{#, fInv[#, f[#, 0.95]]}] &) /@Range[415, 420]]

(Note that there were missing underscores in the function definitions that I've added in.)

I get the same result as you with Mathematica 10.1 (Windows 7):

{218.647, {415, 0.95}},
{326.461, {416, 0.95}},
{222.949, {417, 0.95}},
{339.589, {418, 0.95}},
{80.7641, {419, 0.50088}},
{116.356, {420, 0.50338}}

The good news is that it doesn't take as long to get a wrong answer. (Sorry, couldn't help myself.)

So...Is this a bug or just pushing the function too far? If the latter, should there not be a warning?

POSTED BY: Jim Baldwin

Good, thanks. I will practice pasting code next time. I think this is a bug. I just sent it to Wolfram support and I attached that notebook (btw, I see that one could attach files here too). I noticed the underscores, but I thought I had fixed them by hand.

I am not sure this is pushing the function too far. It is odd, but sometimes I get better results for larger numbers (Mathematica seems to change its mind about the method to use, I recall seeing an example to this effect somewhere in the 'Help'). I have plotted the time it takes to make similar calculations (as a function of n, say) and it does very strange things, not random, but difficult to explain.

I am pretty sure Mathematica can do better. I often work with the regression model functions (LinearModelFit and NonlinearModelFit). These functions can provide tables with statistical summaries of the regression models (p-values for F- and t-statistics, etc.). It is quite challenging (but possible) to reproduce the results in those tables using the CDF and InverseCDF functions, but those Fit functions evaluate quickly and print the tables in a heartbeat. It seems clear that there are statistical packages with more efficient versions of the CDF functions somewhere in the code.

Anyway, thanks for your attention and feedback. Wolfram sent me an automated reply saying they will get back to me within 3 business days. I will update this post if they provide any useful work-around.

OL.

POSTED BY: Otto Linsuain

Thanks again. You are right that it takes less than a tenth of a second. It is also interesting that the message prints when calculating the CDF rather than the inverse (first time I see it that way).

You are also right that my original posting was long.

I still think that if I can improve this with a simple bisection code, then there should be room for further improvement.

Thanks,

OL.

POSTED BY: Otto Linsuain

Hello Jim,

Do you think it could be worthwhile reposting the question as shown below?

Input:

a = InverseCDF[NormalDistribution[0, 1], 0.95];

f[n, cl] := InverseCDF[NoncentralStudentTDistribution[n - 1, Sqrt[n]*a], cl]

fInv[n, k] := CDF[NoncentralStudentTDistribution[n - 1, Sqrt[n]*a], k]

Column[(AbsoluteTiming[{#, fInv[#, f[#, 0.95]]}] &) /@ Range[415, 420]]

Output:

{{{172.416517, {415, 0.95}}},

{{263.043509, {416, 0.95}}},

{{171.102316, {417, 0.95}}},

{{265.421146, {418, 0.95}}},

{{59.659498, {419, 0.50088}}},

{{93.751030, {420, 0.50338}}}}

The last column should be 0.95 all the time, the last two numbers are wrong. No warnings are printed. Too bad I do not know how to paste Mathematica code here more neatly.

Thanks, OL.

POSTED BY: Otto Linsuain
Posted 9 years ago

It takes less than a tenth of a second. I wasn't trying to simplify things. While it was clear from your write-up as to what the general problem was, if I might be so bold, your message rambled a bit and I was just trying to get something more specific to shoot at. You'll just be more likely to get help with concise set of statements (using the "Code Sample" formatting).

POSTED BY: Jim Baldwin
Posted 9 years ago

To make a specific example, does the following capture the issue? (This is with Mathematica v. 10.1)

(* Set values *)
\[Nu] = 4200;
\[Alpha] = 0.95;

(* Get t-value *)
t = N[InverseCDF[NoncentralStudentTDistribution[\[Nu], 0], \[Alpha]], 50]

(* Plug value back into CDF function *)
N[CDF[NoncentralStudentTDistribution[\[Nu], 0], t], 50]

with output

1.64522

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in
Statistics`NoncentralDistributionsDump`z$163348 near {Statistics`NoncentralDistributionsDump`z$163348} = {64.8939}.
NIntegrate obtained 1.746085001809244`15.954589770191005*^6694 and 6.297643634774244`15.954589770191005*^6690
for the integral and error estimates. >>

0.95
POSTED BY: Jim Baldwin

Hello Jim,

Thanks for replying to my post.

Yes, it captures the essence of the issue. But your result is accurate (even if it printed a warning).

You simplified it slightly by making the non-centrality parameter equal to zero, which basically means you are using the central student T distribution (even if you invoke the non-central one). In my case I need the non-centrality parameter to be given by the formula in my original post, which makes it proportional to the Sqrt of the sample size. The central distribution is symmetric about the peak and peaks at zero. The non-central one is not only shifted away from zero, but also becomes asymmetric about the peak.

I notice is that version 10 seems to be behaving similarly (you got a warning and probably a somewhat long computation time). However, you got an accurate result. In my case, I have noticed that the result was inaccurate for sample sizes starting around 419.

Thanks,

OL.

POSTED BY: Otto Linsuain
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract