[ https://issues.apache.org/jira/browse/MATH909?page=com.atlassian.jira.plugin.system.issuetabpanels:commenttabpanel&focusedCommentId=13503793#comment13503793
]
Sébastien Brisard edited comment on MATH909 at 11/26/12 2:19 PM:

Patrick, the arguments of the incomplete beta function are large, so a loss of accuracy is
to be expected. We are working on it (see MATH738), but there is still a lot to do.
Meanwhile, I've used MAXIMA to compute in multiple precision the result you are asking for.
Here is the script I wrote
{code}
kill(all);
load(newton1);
fpprec : 128;
d1 : 10675;
d2 : 501725;
p : 0.025b0;
F(x) := block(
y : d1 * x / (d1 * x + d2),
beta_incomplete(0.5b0 * d1, 0.5b0 * d2, y) / beta(0.5b0 * d1, 0.5b0 * d2)
);
x0 : 0.9733505b0;
x1 : newton(F(x)  p, x, x0, 10**(fpprec+1));
print(x1);
{code}
Basically, {{F(X)}} is the cumulative distribution function, and {{p}} is the target. I then
look for an approximate solution of {{F(X) == p}}, starting from {{x0 = 0.9733505b0}} which
is the value returned by R.
Here is the result I get
{code}
9.730779455235159798387713198169726833852358770410804399004797817813867760559669310857118225960097152689836562586854107391295331b1
{code}
It seems to me that the value returned by CM is more accurate than the one returned by R.
Could you carry out an independent check on that?
NOTA: for some reason, I had to compute the regularized incomplete beta function as the ratio
of the incomplete beta function and the beta function, as the function {{beta_incomplete_regularized
(a, b, z)}} led to errors.
was (Author: celestin):
Patrick, the arguments of the incomplete beta function are large, so a loss of accuracy
is to be expected. We are working on it (see MATH738), but there is still a lot to do.
Meanwhile, I've used MAXIMA to compute in multiple precision the result you are asking for.
Here is the script I wrote
{code}
kill(all);
load(newton1);
fpprec : 128;
d1 : 10675;
d2 : 501725;
p : 0.025b0;
F(x) := block(
y : d1 * x / (d1 * x + d2),
beta_incomplete(0.5b0 * d1, 0.5b0 * d2, y) / beta(0.5b0 * d1, 0.5b0 * d2)
);
x0 : 0.9733505b0;
x1 : newton(F(x)  p, x, x0, 10**(fpprec+1));
print(x1);
{code}
Basically, {{F(X)}} is the cumulative distribution function, and {{p}} is the target. I then
look for an approximate solution of {{F(X) == p}}, starting from {{x0 = 0.9733505b0}} which
is the value returned by R.
Here is the result I get
{code}
9.730779455235159798387713198169726833852358770410804399004797817813867760559669310857118225960097152689836562586854107391295331b1
{code}
It seems to me that the value returned by CM is more accurate than the one returned by R.
Could you carry out an independent check on that?
> FDistribution NoBracketingException in BrentSolver
> 
>
> Key: MATH909
> URL: https://issues.apache.org/jira/browse/MATH909
> Project: Commons Math
> Issue Type: Bug
> Affects Versions: 3.0
> Reporter: Patrick Meyer
> Original Estimate: 24h
> Remaining Estimate: 24h
>
> I get an exception when running the code below. the exception is
> {code}
> function values at endpoints do not have different signs, endpoints: [0, 1.002], values:
[0.025, ∞]
> {code}
> The problematic code:
> {code}
> double df1 = 10675;
> double df2 = 501725;
> FDistribution fDist = new FDistribution(df1, df2);
> System.out.println(fDist.inverseCumulativeProbability(0.025));//NoBracketingException
> {code}
> However, R returns the value 0.9733505. The R code is:
> {code}
> qf(p=.025, df1=10675, df2=501725)
> {code}
> I don't know enough about the FDistribution class to know the solution to the exception,
but I thought I would report it.

This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators
For more information on JIRA, see: http://www.atlassian.com/software/jira
