Stan with Stata, Part IV: Automation

Last time I presented a very basic do file for running Stan from within Stata and although the program worked, there were places where we had to find work-arounds that meant that the process was not fully automatic. In particular, we had to write the script manually, we used Stata’s shell command to run the compiled code and we had to copy and paste the results into Stata. This week I will address those issues.

Executables.txt
In order to tell the wbs programs were to find WinBUGS, OpenBUGS etc., I recommend creating a file called executables.txt in your PERSONAL folder (you can find the location of your PERSONAL folder by typing the sysdir command into Stata). executables.txt is a simple text file, which in my case includes lines such as,
WINBUGS,”C:\Software\WinBUGS14\WinBUGS14.exe”
OPENBUGS,”C:\Software\OpenBUGS\OpenBUGS.exe”
JAGS,”C:\Software\JAGS\JAGS-3.4.0\x64\bin\jags.bat”
Now we need to add a line for Stan. We will not give the name of the executable as that is not the way that Stan works, but rather we will give the location of the cmdstan folder. Notice that Stan needs linux-style slashes / in the path.
STAN,”C:/Software/cmdstan-v2.6.2″

wbsscript
When the book Bayesian analysis with Stata was published wbsscript only handled WinBUGS and OpenBUGS so when I added the facility to run JAGS, I created a new program called wbsscript_v2.ado, which was downloadable through my webpage. The command was given a new name so as not to cause confusion with the version made available by Stata at the time when the book was published. When the new version has been thoroughly tested I plan to update wbsscript with the contents of wbsscript_v2, except that now I also want to add a stan option, so I have created a new program called wbsscript_v3.ado that updates wbsscript_v2.ado and which can handle WinBUGS, OpenBUGS, JAGS and Stan. This program can also be downloaded from my home page although it should be treated as still being under test.

wbsrun & wbscoda
I have created new versions of wbsrun and wbscoda called wbsrun_v3 and wbscoda_v3 that have a Stan option and which are able to run simple Stan programs and read the results into Stata.

Putting these things together we can run Stan much as we have been running WinBUGS, OpenBUGS and JAGS. Here is an updated version of the program that I presented last week.
clear
set obs 20
local TRUE_THETA 0.4
gen y = runiform() < `TRUE_THETA’

wbsmodel thisdofile.do BernModel.stan
/*
data {
    int N;
    int y[N];
}
parameters {
    real theta;
}
model {
    theta ~ beta(1,1);
    for( i in 1:N ) {
       y[i] ~ bernoulli(theta);
    }
}
*/

wbslist_v2 (N <- 20) (vect y , format(%3.0f)) using BernData.txt, replace jags
wbsscript_v3 using BernScript.bat , replace ///
     model(BernModel.stan) data(BernData.txt) stan
wbsrun_v3 using BernScript.bat , stan
type log.txt
wbscoda_v3 using output.csv , stan clear

mcmctrace theta
mcmcstats theta

Options
Many of the command line options for Stan control aspects of the samplers and we cannot really discuss those without first understanding the nature of Stan’s samplers. However, there are simple features such as the run length and the names of the output files that we can easily modify. I have implemented the following using existing options for wbsscript:
burnin=length of burnin (num_warmup in Stan)
updates = length of the run after the burnin  (num_samples in Stan is updates+burnin)
thin=thinning constant
codafile=output file name (defaults to output.csv for Stan. Of course, Stan creates a .csv file that is not in coda format, so ideally we should use a more generic name for the option)
seed= the random number seed. In Stan a negative value draws a seed from the time of day.

So now we can write commands such as,

wbsscript_v3 using BernScript.bat , replace ///
    model(BernModel.stan) data(BernData.txt) ///
    updates(3000) burnin(500) thin(2) ///
    coda(BernOut.csv) stan
wbsrun_v3 using BernScript.bat , stan
wbscoda_v3 using BernOut.csv , stan clear

To make any more progress we will need to know something about the samplers that Stan uses so that we can call Stan in a more intelligent way. I’ll start to explain the Stan sampler in my next posting. Unfortunately (for the flow of the blog) I will be working abroad for the next couple of weeks so the explanation of Hamiltonian Monte Carlo (also known as Hybrid Monte Carlo) will have to wait for my return.

If you want to download any of the programs referred to in this posting they are available from my webpage at https://www2.le.ac.uk/Members/trj