# Random number generation for conjugate distribution of beta distribution

I try to generate random numbers from the conjugate distribution of beta distribution. It is as follows

$$p(α,β∣a,b,d)∝ \frac{e^{-a \alpha} e^{-b \beta}}{(\beta(\alpha,\beta))^d} \:\:\:\:,\:\:\: \alpha>0,\beta>0$$

where $$a>0$$, $$b>0$$ and $$d>0$$. $$\beta(\alpha,\beta)$$ is the Beta function. How can I generate samples from the above distribution? Thank you.

• This answer to another question gives you the proper link. Aug 29, 2020 at 15:46

Here is an excerpt from our book, Introducing Monte Carlo methods with R, indirectly dealing with this case (by importance sampling). The graph of the target shows a smooth and regular shape for the conjugate, meaning a Normal or Student proposal could maybe be of used for accept-reject. An alternative is to use MCMC, eg Gibbs sampling.

Example 3.6. [p.71-75] When considering an observation $$x$$ from a beta $$\mathcal{B}(\alpha,\beta)$$ distribution, $$x\sim \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\,x^{\alpha-1} (1-x)^{\beta-1}\,\mathbb{I}_{[0,1]}(x),$$ there exists a family of conjugate priors on $$(\alpha,\beta)$$ of the form $$\pi(\alpha,\beta)\propto \left\{ \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha) \Gamma(\beta)} \right\}^\lambda\, x_0^{\alpha}y_0^{\beta}\,,$$ where $$\lambda,x_0,y_0$$ are hyperparameters, since the posterior is then equal to $$\pi(\alpha,\beta|x)\propto \left\{ \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha) \Gamma(\beta)} \right\}^{\lambda+1}\, [x x_0]^{\alpha}[(1-x)y_0]^{\beta}\,.$$ This family of distributions is intractable if only because of the difficulty of dealing with gamma functions. Simulating directly from $$\pi(\alpha,\beta|x)$$ is therefore impossible. We thus need to use a substitute distribution $$g(\alpha,\beta)$$, and we can get a preliminary idea by looking at an image representation of $$\pi(\alpha,\beta|x)$$. If we take $$\lambda=1$$, $$x_0=y_0=.5$$, and $$x=.6$$, the R code for the conjugate is

f=function(a,b){
exp(2*(lgamma(a+b)-lgamma(a)-lgamma(b))+a*log(.3)+b*log(.2))}


leading to the following picture of the target: The examination of this figure shows that a normal or a Student's $$t$$ distribution on the pair $$(\alpha,\beta)$$ could be appropriate. Choosing a Student's $$\mathcal{T}(3,\mu,\Sigma)$$ distribution with $$\mu=(50,45)$$ and $$\Sigma=\left( \begin{matrix}220 &190\\ 190 &180\end{matrix}\right)$$ does produce a reasonable fit. The covariance matrix\idxs{covariance matrix} above was obtained by trial-and-error, modifying the entries until the sample fits well enough:

 x=matrix(rt(2*10^4,3),ncol=2)       #T sample
E=matrix(c(220,190,190,180),ncol=2) #Scale matrix
image(aa,bb,post,xlab=expression(alpha),ylab=" ")
y=t(t(chol(E))%*%t(x)+c(50,45))
points(y,cex=.6,pch=19)


If the quantity of interest is the marginal likelihood, as in Bayesian model comparison (Robert, 2001), $$\begin{eqnarray*} m(x) &=& \int_{\mathbb R^2_+} f(x|\alpha,\beta)\,\pi(\alpha,\beta)\,\text{d}\alpha \text{d}\beta \\ &=& \dfrac{\int_{\mathbb R^2_+} \left\{ \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha) \Gamma(\beta)} \right\}^{\lambda+1}\, [x x_0]^{\alpha}[(1-x)y_0]^{\beta} \,\text{d}\alpha \text{d}\beta} {x(1-x)\,\int_{\mathbb R^2_+} \left\{ \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha) \Gamma(\beta)} \right\}^{\lambda}\, x_0^{\alpha} y_0^{\beta} \,\text{d}\alpha \text{d}\beta}\,, \end{eqnarray*}$$ we need to approximate both integrals and the same $$t$$ sample can be used for both since the fit is equally reasonable on the prior surface. This approximation \begin{align}\label{eq:margilike} \hat m(x) = \sum_{i=1}^n &\left\{ \frac{\Gamma(\alpha_i+\beta_i)}{\Gamma(\alpha_i) \Gamma(\beta_i)} \right\}^{\lambda+1}\, [x x_0]^{\alpha_i}[(1-x)y_0]^{\beta_i}\big/g(\alpha_i,\beta_i) \bigg/ \nonumber\\ &x(1-x)\sum_{i=1}^n \left\{ \frac{\Gamma(\alpha_i+\beta_i)}{\Gamma(\alpha_i) \Gamma(\beta_i)} \right\}^{\lambda}\, x_0^{\alpha_i}y_0^{\beta_i}\big/g(\alpha_i,\beta_i)\,, \end{align} where $$(\alpha_i,\beta_i)_{1\le i\le n}$$ are $$n$$ iid realizations from $$g$$, is straightforward to implement in {\tt R}:

 ine=apply(y,1,min)
y=y[ine>0,]
x=x[ine>0,]
normx=sqrt(x[,1]^2+x[,2]^2)
f=function(a) exp(2*(lgamma(a[,1]+a[,2])-lgamma(a[,1])
-lgamma(a[,2]))+a[,1]*log(.3)+a[,2]*log(.2))
h=function(a) exp(1*(lgamma(a[,1]+a[,2])-lgamma(a[,1])
-lgamma(a[,2]))+a[,1]*log(.5)+a[,2]*log(.5))

den=dt(normx,3)

> mean(f(y)/den)/mean(h(y)/den)
 0.1361185


Our approximation of the marginal likelihood, based on those simulations is thus $$0.1361$$. Similarly, the posterior expectations of the parameters $$\alpha$$ and $$\beta$$ are obtained by

> mean(y[,1]*f(y)/den)/mean(f(y)/den)
 94.08314
> mean(y[,2]*f(y)/den)/mean(f(y)/den)
 80.42832


i.e., are approximately equal to $$19.34$$ and $$16.54$$, respectively.