Eu estou procurando uma biblioteca / implementação java que suporta o cálculo da função de distribuição cumulativa inversa para a distribuição beta (aka estimativa de quantis)com razoável precisão.

Claro que eu tenteiapache commons matemática, mas na versão 3 ainda parece haver algumquestões com a precisão. Abaixo, o problema que levou a essa questão é descrito extensivamente.

Suponha que eu queira calcular o intervalo credível de uma distribuição beta com muitos testes. Emapache commons matemática ...

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

A questão é que o percentil 2,5 e a mediana são os mesmos, entretanto, ambos maiores que a média.

Em comparação, oR-pacotebinom 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

e aR-pacoteEstatísticas

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

Para secundar os resultados de R, aqui está o queWolfram Alpha me disse

InverseBetaRegularized [0,025,10007 + 1,161750-10007 + 1] => 0.06070354631 ...InverseBetaRegularized [0.975,10007 + 1,161750-10007 + 1] => 0,06305170794 ...

Nota final sobre os requisitos:

Eu preciso executar muitos desses cálculos. Portanto, qualquer solução não deve demorar mais do que 1s (o que ainda é muito comparado com os 41ms (embora errados) do apache commons math).Estou ciente de que se pode usar R dentro de java. Por razões que não vou elaborar aqui, esta é a última opção se qualquer outra coisa (java puro) falhar.

Atualização 21.08.12

Parece que o problema foi corrigido ou pelo menos melhorado em 3.1-SNAPSHOT de apache-commons-math. Para o usecase acima

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

Atualização 23.02.13

Embora, à primeira vista, essa questão e suas respostas possam ser muito localizadas, acho que isso ilustra muito bem que alguns problemas numéricos não podem ser resolvidos (eficientemente) com uma abordagem do tipo "o que vem primeiro à mente-hacker". Então, espero que continue aberto.