Cómo calcular la función de distribución beta acumulativa inversa en java

Estoy buscando una biblioteca / implementación java que admita el cálculo de la función de distribución acumulativa inversa para la distribución beta (también conocida como estimación de cuantiles)con razonable precisión.

Claro que lo he intentadoapache commons math, pero en la versión 3 todavía parece haber algunosProblemas con la precisión.. A continuación se describe ampliamente el problema que llevó a esta pregunta.

Supongamos que quiero calcular el intervalo creíble de una distribución beta con muchas pruebas. Enapache commons math ...

final int trials = 161750;
final int successes = 10007;
final double alpha = 0.05d;

// the supplied precision is the default precision according to the source code
BetaDistribution betaDist = new BetaDistribution(successes + 1, trials - successes + 1, 1e-9);

System.out.println("2.5 percentile :" + betaDist.inverseCumulativeProbability(alpha / 2d));
System.out.println("mean: " + betaDist.getNumericalMean());
System.out.println("median: " + betaDist.inverseCumulativeProbability(0.5));
System.out.println("97.5 percentile :" + betaDist.inverseCumulativeProbability(1 - alpha / 2d));

que entrega

2.5 percentile :0.062030402074808505
mean: 0.06187249616697166
median: 0.062030258659508855
97.5 percentile :0.06305170793994147

El problema es que el percentil 2.5 y la mediana son iguales mientras que ambos son mayores que la media.

En comparación, laR-paquetebinom entrega

binom.confint(10007+1,161750+2,methods=c("agresti-coull","exact","wilson"))
         method     x      n      mean      lower      upper
1 agresti-coull 10008 161752 0.0618725 0.06070873 0.06305707
2         exact 10008 161752 0.0618725 0.06070317 0.06305756
3        wilson 10008 161752 0.0618725 0.06070877 0.06305703

y elR-paqueteestadísticas

qbeta(c(0.025,0.975),10007+1,161750-10007+1)
[1] 0.06070355 0.06305171

A segundo los resultados de R, esto es lo queWolfram Alpha me dijo

InverseBetaRegularized [0.025,10007 + 1,161750-10007 + 1] => 0.06070354631 ...InverseBetaRegularized [0.975,10007 + 1,161750-10007 + 1] => 0.06305170794 ...

Nota final sobre los requisitos:

Necesito hacer muchos de estos cálculos. Por lo tanto, cualquier solución no debe durar más de 1 s (que aún es mucho en comparación con los 41 ms de las matemáticas de apache commons (aunque mal)).Soy consciente de que uno puede usar R dentro de java. Por razones que no explicaré aquí, esta es la última opción si falla alguna otra cosa (Java puro).

Actualización 21.08.12

Parece que el problema se ha solucionado o al menos mejorado en 3.1-SNAPSHOT de apache-commons-math. Para el caso de uso anterior.

2.5 percentile :0.06070354581340706
mean: 0.06187249616697166
median: 0.06187069085946604
97.5 percentile :0.06305170793994147

Actualización 23.02.13

Si bien a primera vista esta pregunta y sus respuestas pueden estar demasiado localizadas, creo que ilustra muy bien que algunos problemas numéricos no se pueden resolver (de manera eficiente) con un enfoque de lo primero que viene a la mente del hacker. Así que espero que quede abierto.

Respuestas a la pregunta(3)

Su respuesta a la pregunta