Prof Darren Wilkinson Professor of Stochastic Modelling
School of Mathematics & Statistics
Newcastle University

Generating beta random quantities

Let's now look at a simple way of generating Be(a,b) random variates. Note that this is not a very efficient way to do it, particulary for very concentrated distributions, but serves to illustrate the basic method.

We consider only the case a,b>1 so that the density is bounded and has mode at (a-1)/(a+b-2). We can just use the scaled density x^(a-1)*(1-x)^(b-1), and simulate an X on (0,1) and a Y between 0 and the height at the mode. Accept the X if the Y is less than the scaled density at X.

Implementation in R

The logic here is now sufficiently complex that it is worth defining a few auxilliary functions encapsulating the density, its mode, and its maximum value. 
sdbeta<-function (x, alpha, beta) 
{
        (x^(alpha - 1)) * ((1 - x)^(beta - 1))
}

betamode<-function (alpha, beta) 
{
        (alpha - 1)/(alpha + beta - 2)
}

sdbetamode<-function (alpha, beta) 
{
        sdbeta(betamode(alpha, beta), alpha, beta)
}
Once we've defined these, the function to simulate a single observation can be written as 
myrbeta1<-function (alpha, beta) 
{
        repeat {
                x <- runif(1, 0, 1)
                y <- runif(1, 0, sdbetamode(alpha, beta))
                fx <- sdbeta(x, alpha, beta)
                if (y < fx) 
                        return(x)
        }
}
and then vectors can be generated in the usual way.

The R source code is available here: rbeta.r

Implementation in C

We can use exactly the same logic in C: 
float sdbeta(float x, float alpha, float beta)
{
  return(pow(x,alpha-1)*pow(1-x,beta-1));
}

float betamode(float alpha, float beta)
{
  return((alpha-1)/(alpha+beta-2));
}

float sdbetamode(float alpha, float beta)
{
  return(sdbeta(betamode(alpha,beta),alpha,beta));
}

float genbeta(float alpha, float beta)
{
  float x,y,fx;
  while (1)
    {
      x=urand(0,1);
      y=urand(0,sdbetamode(alpha,beta));
      fx=sdbeta(x,alpha,beta);
      if (y<fx)
        break;
    }
  return(x);
}
and then call these from a main procedure in the usual way.

The full C source for this is here: beta.c

Other implementations

Those who are interested may also wish to look at corresponding code for LISP-STAT: beta.lsp, Python: beta.py and Sather: beta.sa.

Darren J Wilkinson Book: Stochastic Modelling for Systems Biology  
darren.wilkinson@ncl.ac.uk
http://www.staff.ncl.ac.uk/d.j.wilkinson/