Generate beta random number without statistics toolbox
5 views (last 30 days)
Show older comments
Rémi BERNARD
on 18 Jun 2019
Commented: John D'Errico
on 26 Jun 2019
Hello everyone,
I am currently working on a project involving the PDF of beta distribution, and generation of random numbers with according distribution. My script is supposed to be converted into an executable so I am limited in use of toolboxes and some commands.
Does someone know how to generate such numbers (a few millions of them, so manually is excluded) without the Statistics toolbox?
Many thanks!
0 Comments
Accepted Answer
John D'Errico
on 18 Jun 2019
Edited: John D'Errico
on 18 Jun 2019
As a followup, I note that the answer by dpb seems to be incorrect.
There are several ways to generate a beta random variable. Simplest is the direct use of betaincinv. (The hardest part of that is spelling betaincinv.) What I do not know is if you can use that utility, in an executable. It should be possible.
For example, let me pick alpha = 2, beta = 3, with a sample size of 1e6.
alpha = 2;
beta = 3;
N = 1e6;
U = rand(N,1);
tic,Z = betaincinv(U,alpha,beta);toc
Elapsed time is 0.432232 seconds.
histogram(Z,100,'norm','pdf');
hold on
H = fplot(@(x) betapdf(x,alpha,beta),[0,1]);
H.Color = 'r';
That is clearly reasonable. So if you can use betaincinv, then you are done. As I show in my comment on the answer by dpb, you can also use the ratio from a couple of gamma random variables. But betaincinv seems easiest. It seems reasonably fast given what it does, at less than .5 seconds for a million such events.
And if you cannot use betaincinv, then gammaincinv would also fail. So as long as you can use betaincinv, use that.
4 Comments
John D'Errico
on 26 Jun 2019
In fact, it is the same distribution as what betarnd would yield (unless alpha and beta are very near zero. That is itself an issue that needs discussion, later.)
Why do I use rand in there, and then call betaincinv? That gets into the theory of how one generates random variables besed on other distributions. (There are books to found on the subject.) There are multiple ways to do so, but as long as the inverse function for the distribution CDF is available and computable, then the method I showed is correct. And the nice thing is, the beta CDF is intimately related to the incomplete beta function itself, and betaincinv is exactly what is needed to compute the inverse of the beta distribution CDF.
In the case of very small alpha/beta, you will have a problem with any method explained here. So unless that is an issue, I would just use the betaincinv method, as simple.
An alternative, of computing the ratio X/(X+Y), where X and Y are GAMMA random variables, is an option. I recall dicussing that somewhere in a comment. TWO calls to gammaincinv will then be required, but they may be more efficient than one call to betaincinv. That would take some testing to resolve.
Finally, is the case where alpha and beta are both fairly small numbers. In that event, the resulting beta distribution becomes fairly nasty, and strongly bimodal. So one could then be forced to use a binomial distribution, as a special case. Only you might know if that is a case to worry about.
More Answers (1)
dpb
on 18 Jun 2019
Given Ui are uniform random variates
Y1=U1/α1 and Y2=U1/β1,
such that Y1+Y2<=1,
then
X=Y1/(Y1+Y2)
follows the beta distribution with parameters α and β
5 Comments
John D'Errico
on 19 Jun 2019
As I said, what dpb told you to do is not technically correct. In fact, it does not generate the desired distribution, as I proved that it cannot. But you will need to make that decision for yourself.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!