4
$\begingroup$

I know that we can approximate Beta distribution to Normal distribution when the values of alpha and beta are large numbers.

In my problem alpha lies between 1 and 10, beta is always greater than 1000. Is there any distribution (for such a combination of alpha and beta) that we can approximate it to, with a simpler closed form so that it will be easy to calculate CDF and InverseCDF?

(I have tried Kumaraswamy distribution but did not get good results due to the large values of beta.)

EDIT "easy"- in terms of computation time in JAVA/PHP.

(I need an approximate distribution since it takes lot of time to compute the InverseCDF of Beta distribution(0.05 ms in java is also very slow for my needs).

For alpha greater than 10 i approximated Beta distribution to Normal distribution and computed InverseCDF like this which included just math. So i did not use any other library functions like Normal distribution or Beta distribution, which reduced the processing time drastically.)

Improve this question
$\endgroup$
12
  • 2
    $\begingroup$ Why do you need an approximation? Why not just use the exact beta distribution? $\endgroup$ – Greg Snow Oct 3 '13 at 17:19
  • $\begingroup$ The critical thing in your case is the very large difference between $\alpha$ and $\beta$. Have you tried to graph the theoretical pdf-cdf?. Stick with the beta distribution itself. $\endgroup$ – Alecos Papadopoulos Oct 3 '13 at 18:13
  • 2
    $\begingroup$ @Alecos Notice that the PDF of the Beta$(\alpha,\beta)$ distribution, $x^\alpha(1-x)^{\beta-1}\frac{dx}{x},$ can be rescaled by $1/\beta$ to yield $(\frac{x}{\beta})^\alpha(1-\frac{x}{\beta})^{\beta-1}\frac{dx}{x}$ and that for large $\beta$, $(1-\frac{x}{\beta})^{\beta-1}\approx\exp(-x)$, giving the PDF for a scaled Gamma$(\alpha)$ distribution. By matching expectations you get a remarkably good approximation; for $\alpha\le 10$, the absolute difference in CDFs is comparable to $1/\beta$. So maybe there is some merit in an approximation if one finds Gamma integrals easy to compute. $\endgroup$ – whuber♦ Oct 3 '13 at 18:18
  • 3
    $\begingroup$ My point, @Alecos, is that your expectation of huge errors is not borne out: the approximation is excellent in all respects, whether measured in terms of the PDFs, the CDFs, or lower moments. The message is that a Beta with extreme values like these is practically a Gamma, in much the same way that a Binomial with extreme values of its parameter is practically a Poisson. In principle, then, using the Gamma approximation should have uses similar to applying Poisson models in problems that, sensu strictu, are really Binomial. $\endgroup$ – whuber♦ Oct 3 '13 at 18:31
  • 3
    $\begingroup$ It's just that your knowledge and understanding of statistics are outstanding.... $\endgroup$ – Nick Cox Oct 4 '13 at 17:24
6
$\begingroup$

Johnson and Kotz point out that

  1. The $B$ (Beta) and $F$ distributions are closely related: When $k$ has a $B(\nu_1, \nu_2)$ distribution, then $f$ has an $F(2\nu_1, 2\nu_2)$ distribution where

    $$f = \frac{k \nu_2}{(1-k)\nu_1}.$$

  2. The Wilson-Hilferty approximation to $\chi^2$ (or, equivalently, a $\Gamma$ distribution)--which has "remarkable accuracy"--can be applied to $F$, which itself is a ratio of $\chi^2$ distributions.

The result is that

$$z = \frac{\sqrt[3]{f} \left(1-\frac{2}{9 \nu_2}\right)+\left(\frac{2}{9 \nu_1}-1 \right)}{\sqrt{\frac{2 f^{2/3}}{9 \nu_2}+\frac{2}{9 \nu_1}}}$$

has approximately a standard Normal distribution. This provides relatively fast calculations of the CDF of either the $B$ or the $F$ distribution and is not difficult to invert.

This approximation works fairly well for $\nu_2 \ge 1000$ and $\nu_1=1$ and progressively better as $\nu_1$ increases. Here is a set of plots of the difference between the approximate CDF and the correct CDF for $\nu_2=1000$ and $\nu_1 = 1,2,3,6,10$, shown with identical scales on their axes:

CDF error plots

The error is extremely low in the right tail, where interest often focuses, but even in the left tail it rarely exceeds $0.002$ once $\nu_1 \ge 2$. As $\nu_2$ increases, the errors remain about the same.

Given the level of accuracy of this approximation, it makes little sense to compute the standard Normal CDF (or its inverse) to double-precision accuracy. This will allow for even faster computation.

One approach is to create a table of values and interpolate within it. Because the Normal CDF is so curved near the tails, it would seem a largish table might be needed for high accuracy. However, by considering Mills Ratio one might be led instead to create a (very) small table of pairs $(z, \sqrt{-2\log(\Phi(z))}$ for $z$ in the lower tail, say for $-8 \le z \le 0$, where $\Phi$ is the standard Normal CDF. Cubic spline interpolation does a great job even when these pairs are evaluated only for integral $z$ (nine values total in the entire table). To further simplify the programming (if a cubic spline interpolation routine is not readily available) one could implement a quadratic spline. Even better--and barely more difficult--is to form the weighted average of a linear interpolator (which tends to underestimate the CDF) and a quadratic interpolator (which overestimates the CDF). A weight of $0.238$ applied to the linear interpolator works well:

Comparison with normal CDF

This plot shows the differences between three interplators and the standard Normal CDF: the red one is for the linear interpolator, the gold for the quadratic interpolator, and the blue (in between) for their weighted average.

(The same technique of averaging in a linear interpolator also decreases the error by an order of magnitude when using a cubic spline.)

The weighted average is another quadratic interpolator but it is no longer a spline (it fails to have continuous first derivatives at the negative integers)--but so what? In effect, it computes the CDF with absolute accuracy better than $0.00005$ for $z\le -1$ and, by symmetry, for $z\ge 1$. Quadratic or cubic interpolation between a small number of pairs of $(z, \Phi(z))$ should work fine for $-1 \lt z \lt 1$ if such calculations are needed. (Typically only very low accuracy, say to about $0.01$, is needed in this range anyway.) The upshot is that the CDF (or, using similar methods, its inverse) can be computed to sufficient accuracy using about a dozen simple arithmetic operations and one exponentiation, at a cost of storing around a dozen precomputed values: even on an interpreted platform like Java that should go very quickly, taking a small fraction of a microsecond.

Examples

Mathematica code to produce the first figure follows.

fDist[f_, {ν1_, ν2_}] := CDF[NormalDistribution[]][((1-2/(9 ν2)) f^(1/3) - (1-2/(9 ν1))) /
    Sqrt[2 f^(2/3)/(9 ν2) + 2/(9 ν1)]];
f[k_, {ν1_, ν2_}] := (k ν2)/((1 - k) ν1);
With[{ν1 = 1, ν2 = 1000},
 Plot[fDist[f[k,{ν1,ν2}], {2 ν1,2 ν2}] - CDF[FRatioDistribution[2 ν1,2 ν2]][f[k,{ν1,ν2}]], 
    {k, 0, 10 ν1/(ν1 + ν2)}, PlotRange -> Full]
 ]

Code to perform the interpolation for $|z| \ge 1$ and plot its error:

d = Table[{z, Sqrt[-2 Log[CDF[NormalDistribution[], z]]]}, {z, -8, -0, 1}];
f1 = Interpolation[d, InterpolationOrder -> 1];
f2 = Interpolation[d, InterpolationOrder -> 2];
f[z_] := 0.238 f1[z] + (1 - 0.238) f2[z];
Plot[{Exp[-f[z]^2/2 ] - CDF[NormalDistribution[], z]}, {z, -4, -1}, PlotRange -> Full]

References

Norman L. Johnson and Samuel Kotz, Continuous Univariate Distributions - 2. J. Wiley & Sons, New York, 1970.

Improve this answer
$\endgroup$
3
  • $\begingroup$ +1, great answer. At first, I didn't see how you have arrived at the formula for the approximation of the $F$ distribution using the Wilson-Hilferty transformations for $\chi^{2}$ random variables. I have found the explanations at the beginning of this paper quite useful. I wonder if you have used the same reasoning (using Fiellers' result) or a completely different derivation? $\endgroup$ – COOLSerdash Oct 5 '13 at 16:42
  • $\begingroup$ @whuber Thank you for the answer :) But i din't clearly understand " The result is z= '' has approximately a standard Normal distribution". Can you please explain that point little bit. (It might be a very trivial thing but i never studied or did anything related to Chi-square or F distributions) $\endgroup$ – Karthik Oct 8 '13 at 9:23
  • $\begingroup$ It means that we may approximate the chances of any event of the form $z\in A$ by the integral $\int_A\exp(-x^2/2)dx/\sqrt{2\pi}.$ This reduces the problem to computing such integrals. The additivity property of integrals implies it suffices to compute probabilities of the form $F(a)=\Pr(z\le a)$ for arbitrary numbers $a$. $F$ is called the cumulative distribution function (CDF). It is often written $\Phi$ for the standard Normal distribution, whence $\Phi(a)=\int_{-\infty}^a\exp(-x^2/2)dx/\sqrt{2\pi}.$ My analysis shows how to approximate these integrals to speed up the calculation. $\endgroup$ – whuber♦ Oct 8 '13 at 14:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.