[
https://issues.apache.org/jira/browse/MATH-909?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13503793#comment-13503793
]
Sébastien Brisard edited comment on MATH-909 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 MATH-738), 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.730779455235159798387713198169726833852358770410804399004797817813867760559669310857118225960097152689836562586854107391295331b-1
{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 MATH-738), 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.730779455235159798387713198169726833852358770410804399004797817813867760559669310857118225960097152689836562586854107391295331b-1
{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: MATH-909
> URL: https://issues.apache.org/jira/browse/MATH-909
> 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
