MCMC algorithms can be slow, so it is often necessary to pay particular attention to the efficiency of one’s code and usually this means programming in Mata. For this reason, the slice, griddy, ARS and ARMS samplers that are described in ‘Bayesian Analysis with Stata’ were programmed using Mata even though they can be called from either Mata or Stata. For convenience, those four samplers are collected together with other Mata code referred to in the book, in a library that I called libmcmc (in Stata a library name needs to begin with lib).
In future postings I will introduce Mata code for fitting some further Bayesian models, so I thought that I would create a new library that I will call libbayes. This library will contain all of the Mata functions that I refer to in this blog. I will start today by showing how the library is created and then I will add some functions for simulating random values from a few standard distributions.
In order to create the library, we need a do file that follows the pattern:
set matastrict off
mata:
mata clear
mata mlib create libbayes , dir(PERSONAL) replace
void prog1(…)
{ ….
}
mata mlib add libbayes prog1()
void prog2(…)
{ ….
}
mata mlib add libbayes prog2()
….
mata mlib index
end
When this do file is run, the library is created in the PERSONAL folder, which I have chosen because that directory is searched automatically by Stata when it looks for Mata libraries (type the command -sysdir- to find the location of this folder). Each program is added in turn and then Stata is told to update its index of Mata libraries. Updating the index is only needed the first time that a new library is created because the index is automatically updated at the start of each Stata session, but it does no harm if the index is updated unnecessarily. I tend to be rather lazy with my programming so I have switched matastrict off. Had I switched it on, then I would have had to define the type of every structure that I use in the functions.
As my first project I plan to write a function for fitting a finite mixture of multivariate normal distributions using conjugate priors and a Gibbs sampler. This code can be used in many different ways, for example for density smoothing, clustering or even fitting smooth lines to a scatter plot, it also serves as a convenient stepping stone to introducing non-parametric Bayesian analysis, which can be though of as an extension of finite mixture modelling in which the number of components in the mixture becomes infinite.
To create the Gibbs sampler for fitting a mixture of multivariate normal distributions, we will need functions that simulate random values from the Wishart, Multivariate normal, Dirichlet and Categorical distributions, so today I will add those functions to my library.
set matastrict off
mata:
mata clear
mata mlib create libbayes , dir(PERSONAL) replace
/*--------------------------------------------
* Single categorical observation
* P = probabilities of each cell (sum to 1)
* returns selected cell number
*--------------------------------------------*/
real scalar rCategorical(real vector P)
{
real scalar k,u,h,sumP
k=length(P)
u=runiform(1,1)
h=1
sumP=P[1]
while( u > sumP & h < k ) {
h++
sumP=sumP+P[h]
}
return(h)
}
mata mlib add libbayes rCategorical()
/*--------------------------------------------
* Single Dirichlet distribution - parameter alpha
*--------------------------------------------*/
real colvector rDirichlet(real vector alpha)
{
real scalar h,sP
real colvector P
h=length(alpha)
P=J(h,1,0)
sP=0
for(i=1;i<=h;i++) {
P[i]=rgamma(1,1,alpha[i],1)
sP=sP+P[i]
}
return(P/sP)
}
mata mlib add libbayes rDirichlet()
/*--------------------------------------------
* Single random multivariate normal variable* M = mean
* L = cholesky(V), V=variance matrix
*--------------------------------------------*/
real rowvector rMNormal(real rowvector M,real matrix L)
{
return(M+L*rnormal(rows(L),1,0,1))
}
mata mlib add libbayes rMNormal()
/*--------------------------------------------
* Single Wishart matrix
* L = cholesky(invsym(R)), where R = k*variance matrix
* k = degree of freedom
*--------------------------------------------*/
real matrix rWishart(real matrix L,real scalar k )
{
p = rows(L)
B = J(p,p,0)
for(i=2;i<=p;i++) {
for(j=1;j<i;j++) {
B[i,j] = rnormal(1,1,0,1)
}
}
for(i=1;i<=p;i++) B[i,i] = sqrt(rchi2(1,1,k-i+1))
return(L*B*B'*L')
}
mata mlib add libbayes rWishart()
mata mlib index
end
The multivariate normal simulator takes the cholesky decomposition of the variance matrix, V, as its input, so it could be called as,
Y = rMNormal(Mu,cholesky(V))
The advantage of calculating the cholesky decomposition outside of the function is that if you want to simulate many values from distributions with the same variance matrix, then you can calculate the cholesky decomposition once and then call rMNormal() as many times as you wish.
A similar reasoning leads us to use cholesky(invsym(R)) as the input to the Wishart function. What is more, anyone who likes to parameterize in terms of S=invsym(R) can use the same function by calling it with the argument cholesky(S).
The Dirichlet distribution is a generalization of the Beta distribution and provides a simple prior for a set of probabilities that sum to one. Drawing random probabilities from a Dirichlet distribution is just a matter of creating a set of appropriately chosen gamma variables and then normalizing them to sum to one. This algorithm together with lots of other information on the Dirichlet distribution can be found on the Wikipedia page (http://en.wikipedia.org/wiki/Dirichlet_distribution )
Testing the functions
The functions in our library will be called by numerous other programs so it is vitally important that we test them to ensure, first, that they are correct and second, that they are robust to misuse. Robustness usually involves checking the inputs to ensure that dimensions match and constraints are fulfilled, for instance, the scalar parameter of a Wishart distribution should not be less than the dimension of the matrix. I have deliberately chosen not to include any such checks in these functions because in a Bayesian analysis they might be called hundreds of thousands of times and speed is important. I leave it to the calling routines to ensure that the parameters are appropriate.
To check the accuracy of the functions we could write a series of Mata programs that call them but I prefer the interactivity of Stata and so I wrote four short Stata programs that call these Mata functions and I used those for testing. Here is the code that I used to test the rCategorical() function.
program rCat syntax newvarlist(max=1) , N(integer) P(namelist) if `n' > _N qui set obs `n' qui gen `varlist' = . mata: P = st_matrix("`p'") forvalues i=1/`n' { mata: H = rCategorical(P) mata: st_numscalar("r(h)",H) qui replace `varlist' = r(h) in `i' } end
matrix PR = (0.2,0.4,0.1,0.3) rCat x , n(1000) p(PR) tabu x
The Stata program rCat calls rCategorical() to generate n values from a single catergorical distibution with probabilities p (which must sum to one, although this is not checked). Inside the Stata program the probabilities are transferred to Mata using the st_matrix() function and then rCategorical() is called repeatedly and the results are returned as locals. After the Stata function has been called, the results are tabulated so that we can see that the proportions agree with the probabilities supplied to the function. This approach is not as efficient as creating a calling function in Mata but it is simple to write and simple to use.
Downloading the code
The do file for creating this library can be downloaded from makeLibbayesVersion1. The software for this blog will not allow me to upload a do file directly, so the link is to a pdf; however, it is simple to copy and paste the code from the pdf into the do file editor. The test programs can be downloaded from testprograms.
Recent Comments