commons-issues mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Sébastien Brisard (JIRA) <j...@apache.org>
Subject [jira] [Comment Edited] (MATH-909) FDistribution NoBracketingException in BrentSolver
Date Mon, 26 Nov 2012 14:20:58 GMT

    [ 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

Mime
View raw message