help amcmc

Title

amcmc -- Mata Functions and Structures for Adaptive Markov Chain Monte Carlo sampling from distributions

Introduction

amcmc refers to a collection of Mata tools for performing adaptive Markov chain Monte Carlo (MCMC) sampling from distributions. Algorithms are
described in Baker (2013) and rely on Markov chain Monte Carlo sampling using a continually adapted multivariate normal proposal
distribution. Using the commands and functions, a user can perform adaptive MCMC using either the function amcmc()
or by setting up an adaptive MCMC
problem using the suite of commands amcmc_*(). The function amcmc() is described below under the heading
of the function amcmc(). Setting up a structured object for performing amcmc_*()
is discussed under the heading using amcmc_*() commands.

The function amcmc()

Syntax

real matrix
amcmc(alginfo,
lnf,
xinit,
Vinit,
draws,
burn,
damper,
aopt,
arate,
vals,
lam,
blocks [,
M,
noisy])

where inputs are

alginfo: string rowvector

lnf: pointer (real scalar function) scalar f

xinit: real colvector

Vinit: real matrixvector

draws: real scalar

burn: real scalar

delta: real scalar

aopt: real scalar

arate: transmorphic

vals: transmorphic

lam: transmorphic

blocks: real matrix

M: transmorphic

noisy: string scalar

Description

amcmc() performs adaptive MCMC sampling using a multivariate normal proposal distribution. MCMC methods work through
acceptance-rejection sampling; an MCMC algorithm with adaptation is continually
tuned as the algorithm runs so as to achieve a targeted acceptance rate. The degree of tuning
recedes as the algorithm proceeds so as
to achieve an eventually stationary proposal distribution. For a full description of the Mata implementation, see Baker (2013),
which follows the more detailed descriptions of adaptive MCMC algorithms in Andrieu and Thoms (2008).
The draws from the target distribution are returned as a matrix, with each row representing a draw.

Discussion of alginfo, the first argument, is postponed until other arguments have been discussed.

The (pointer) function lnf specifies the target distribution from which the user wishes to draw. The function describing the target
distribution must be written in log form, so that the (typically, but not always scalar) log-value of the distribution is returned by the function.

The initial values used in drawing are specified in the argument xinit, and an initial covariance matrix for the
proposal distribution is specified in Vinit. draws tells
the algorithm how many draws to perform, while burn instructs amcmc() to drop the first burn draws as a burn-in period;
accordingly,
amcmc() returns only the last draws-burn draws from the distribution lnf.

delta is an adjustment parameter telling the algorithm how agressively or conservatively to adapt the proposal distribution to achieve the
acceptance rate specified by the user in aopt. delta should lie between zero and one, with values closer to zero corresponding
with slower adaptation, and values closer to one corresponding to more rapid adaptation to the proposal history.
aopt is the acceptance rate desired by the user; typically lying in the range .2 to .44, as
it has been shown that optimal acceptance rates for univariate problems are around .44, while for large-dimensional problems
optimal rates are around .234. See Andrieu and Thoms (2008) for further discussion.

arate, vals, and lam are arguments that can be initialized as missing or as anything else by the user, as they are overwritten.
arate is overwritten with acceptance rates of the algorithm.
vals is overwritten with the values of the target distribution corresponding with each draw.
lam is a set of scaling parameters which are tuned as the algorithm proceeds. These parameters
scale the proposal covariance matrix to direct the algorithm towards the desired acceptance rate, with
aggressiveness captured by delta. In
global, all-at-once sampling, lam is returned as a scalar,
while if metropolis-within-gibbs sampling is specified (so that each component of the
target distribution is drawn sequentially), lam returns a vector
of lambda values equal in dimension to the target distribution. Finally, if block sampling is specified, lam returns a vector of lambda values
equal in dimension to the number of blocks. For further description of how different types are specified, see the description of alginfo below
and the examples. In block sampling or in metropolis-within-gibbs sampling, the dimension of arate matches that of sampler.

blocks is a matrix of zeros and ones which describes how the algorithm is to proceed if the user wishes to draw from the target distribution
not all-at-once but in a sequence of (Gibbs) steps. Values that are to be drawn together are marked by rows of blocks containing ones and zeros
elsewhere.

M can be used to relay additional information to the algorithm. In cases in which the user has assembled a
model statement using moptimize() or optimize(), this model can be passed to amcmc() using M. The idea is to
save the user the need to respecify things such as missing values, arguments of the function embedded in the model, etc. when switching or using amcmc().
In cases in which the user does not have a model, but has a function requiring additional arguments, M can be a pointer holding additional arguments
of the target distribution; up to ten additional arguments can be passed to amcmc() in this fashion. For example, if the target distribution is characterized by
some function lnf(x,Y,Z) where x are to be drawn, but Y and
Z are also arguments, the user would define a pointer M containing Y and Z.
M is optional and does not require specification.

noisy is a string scalar that can be specified as "noisy" or as something else. If noisy="noisy", the algorithm produces
feedback - each time the target distribution is evaluated, it produces a . as output, while after 50 calls, it produces the value of the target
distribution after the last call of lnf.

Last but not least, the first argument of the function is alginfo, which is a string scalar specifying the drawing
scheme desired by the user, what sort of target distribution
evaluator the user has passed to the function when the target is part of a previously specified model statement,
and perhaps some other things about how the algorithm is to be executed. While the examples present
more details, a user may assemble a string row vector composed of one entry from each of the following:

Sampling information mwg, global, block
Model definition moptimize, optimize, standalone
Model evaluator type d*, q*, e*, g*, v*
Other information fast

Thus, if the user wishes to perform Metropolis-within-Gibbs sampling (each component of lnf sampled alone and in order), and has previously
modeled the target distribution as part of a structure using moptimize() and a type d0 evaluator,
the user would specify
(in any order) alginfo="moptimize","d0","mwg". Note that each component of alginfo should be a separate
string entry in a string row vector, so alginfo="moptimize,d0,mwg" will not work. The final option,
"fast", is somewhat experimental
and should be used with caution. In large problems, global, all-at-once samplers require Cholesky decomposition of the proposal covariance matrix.
This can be slow and time consuming. The option "fast" avoids Cholesky decomposition by working with a diagonal proposal covariance matrix that
is continually adjusted as the algorithm proceeds. "fast" sometimes works if the target distribution has many independent random variates and
if the proposal distribution is close to the target distribution. See the examples for other possibilities.

If the user wishes to use a block sampler (which is somewhere between a one-at-a-time, Metropolis-within-Gibbs sampler,
the user must also submit a block matrix to communicate to amcmc() which values are to be drawn together.

Further Options

While the function amcmc() is designed to perform adaptive MCMC, adaptation can be turned off by the user. This feature is often useful if
a previous set of draws from a target distribution has already been performed and the distribution is well-tuned. In this case, the user can set
damper equal to missing, i.e., damper=., and no adaptation will occur. In this case, the user must also specify a set of conformable values
for lam.

Examples

Example 1: Sampling from a univariate normal mixture - with probability 1/2 the mean is -1 or 1
with standard deviation 1 in each case. While this is probably not the most efficient way to sample
from a mixture distribution, it serves to illustrate basic ideas. Initial values for
drawing are set to 0, with initial variance matrix for proposals set at 1.
1000 draws are taken with the first 100 discarded. A value of delta=1/2 is a
fairly conservative choice. A type of sampler has not been specified, which means a global sampler will be used as a default:

:real scalar mixnorm(X)
>{
>val=1/2*normalden(x,1,1)+1/2*normalden(x,-1,1)
>return(ln(val))
>}

:alginfo="standalone"

:X=amcmc(alginfo,&mixnorm(),0,1,1000,100,1/2,.44,arate=.,vals=.,lam=.,.)

Example 2: Sampling from a bivariate normal mixture with dimension 2, where with probability p
the mean is m, and with probability p the mean is -m. p, m and Sig - the covariance matrix of the distribution -
are passed to amcmc() as additional arguments of the function. A vector of zeros and an identity matrix are used as the starting values for the proposal
distribution:

real scalar mixnorm2(x,p,m1,m2,Sig)
>{
>dSig=1/sqrt(det(Sig))
>Siginv=invsym(Sig)
>val1=1/2*dSig*exp(-(x-m1)*Siginv*(x-m1)')
>val2=1/2*dSig*exp(-(x-m2)*Siginv*(x-m2)')
>return(ln(p*val1+(1-p)*val2))
>}

:p=1/2

:m1=1,1

:m2=-1,-1

:Sig=I(2):+.1

:Args=J(4,1,NULL)

:Args[1]=&p

:Args[2]=&m1

:Args[3]=&m2

:Args[4]=&Sig

:alginfo="standalone","global"

:X=amcmc(alginfo,&mixnorm2(),J(1,2,0),I(2),100000,10000,1,.34,arate=.,vals=.,lam=.,.,Args)

Example 3: "Estimating" a negative-binomial count model by simulation. The premise behind the example
is that one can view a likelihood function as a distribution for
parameters conditional on data. The idea is superficially Bayesian; the model is isomorphic to a typical Bayesian
analysis with uninformative prior distributions for parameters, but can be applied more generally; see Chernozukov and Hong (2003).
Sampling from the parameters' distribution is a sometimes useful method for analyzing models applied to
small samples, as results do not depend upon asymptotics as
they do in the typical Bayesian analysis. While in a more complete analysis, one might "thin" the resulting draws out,
check for convergence, etc. this is not done in the example, which is really used to illustrate how a model statement is
passed to amcmc(). The results can be contrasted with those obtained from
estimation via application of nbreg.

. clear all

. use http://www.stata-press.com/data/lf2/couart2, clear
(Academic Biochemists / S Long)

. set seed 5150

. mata:

mata (type end to exit)

:void nb_d0(M,todo,b,lnf,g,H)
>{
>y=moptimize_util_depvar(M,1)
>mu=exp(moptimize_util_xb(M,b,1))
>a=exp(moptimize_util_xb(M,b,2))
>lnfi=lngamma(y:+1/a):-lngamma(1/a):-
> lnfactorial(y):-(y:+1/a):*ln(1:+a*mu):+
> y*ln(a):+y:*ln(mu)
>lnf=colsum(lnfi)
>}

:M=moptimize_init()

:moptimize_init_evaluator(M,&nb_d0())

:moptimize_init_evaluatortype(M,"d0")

:moptimize_init_depvar(M,1,"art")

:moptimize_init_eq_indepvars(M,1,"fem mar kid5 phd ment")

:moptimize_init_eq_indepvars(M,2,"")

:moptimize_evaluate(M)

:alginfo="global","moptimize","d0"

:X=amcmc(alginfo,&nb_d0(),J(1,7,0),I(7),20000,10000,3/4,.234,arate=.,vals=.,lam=.,.,M)

:mean(X)'
1
+
+
1 | -.2200206506 |
2 | .148970965 |
3 | -.1774053762 |
4 | .0173632397 |
5 | .0290353224 |
6 | .253103149 |
7 | -.7998843153 |
+
+

Using amcmc_*() commands

Syntax

Syntax is discussed in a three steps, which include:

Step 1: Initialization
Step 2: Problem definition
Step 3: Obtaining results

Step 1: Initialization

To initialize a problem, the user first issues an initialize command:

A = amcmc_init()

amcmc_init() sets up an adaptive MCMC problem with the following defaults:
the number of draws is set to one, the burn period is set to zero, and the optimal acceptance rate is set to .234.

Step 2: Problem definition

Problem definition fills in the arguments described under the function amcmc(), so it might
first help the user to peruse that material, in addition to having a look at Baker (2013) and/or Andrieu and Thoms (2008). The user has to
relay information about the target distribution, initial values, how many draws are desired, any additional information about the function or
how drawing should go, and how adaptation of the proposal should occur.

amcmc_lnf(A, pointer (real scalar function) scalar lnf)

Use: sets the target distribution as the function lnf. Note that lnf must be in log form.

amcmc_args(A, pointer matrix Args)

Use: sets any additional arguments of the function lnf. Note that this option can only be used with a stand-alone target
distribution. That is,
if the user is sampling from a problem with details constructed using moptimize() or optimize(), the user should
pass any arguments of the function lnf through those routines.

amcmc_xinit(A, real rowvector xinit)

Use: sets initial values for the proposal distribution.

amcmc_Vinit(A, real matrix Vinit)

Use: sets initial values of the proposal covariance matrix. If the matrix submitted is not positive definite, by default a conformable identity
matrix is used.

amcmc_aopt(A, real scalar aopt)

Use: sets desired acceptance rate, typically in the range .234 to .45 or so.

amcmc_damper(A, real scalar delta)

Use: Specifies the parameter determining how tuning of the algorithm is to occur; the value should be between zero and one.
Values close to zero mean less aggressive tuning of the proposal distribution, while values closer to one mean more aggressive tuning. One can
also specify a missing value for damper here; i.e., amcmc_damper(A,.). In this case, no adaptation of the proposal distribution occurs,
and the user must specify scaling parameters using the command amcmc_lambda(A, real rowvector lambda).

amcmc_lambda(A, real rowvector lambda)

Use: Specifies the scaling parameters for the proposal covariance matrix. That is, when a draw is performed using covariance matrix W,
draws use lambda*W to make the draws. This option should only be set if damper (discussed above under the heading
amcmc_damper(A,) is set to missing so that no adaptation of the proposal distribution is to occur.

amcmc_burn(A, real scalar burn)

Use: sets the length of the burn-in period, for which information about draws is discarded.

amcmc_draws(A, real scalar draws)

Use: Number of draws to be performed.

amcmc_noisy(A, string scalar noisy)

Use: If noisy="noisy", a dot is produced each time lnf is called; every 50 calls the function value at the last draw is also
produced.

amcmc_model(A, transmorphic M)

Use: Can be used to append the drawing problem with a previously assembled model statement
formulated using either moptimize() or optimize().

amcmc_blocks(A, real matrix blocks)

Use: In conjunction with block samplers. The matrix blocks contains information about the sequence of draws.

amcmc_alginfo(A, string rowvector alginfo)

Use: A sequence of strings that contain information about how drawing is to proceed, what type of sampler to use, if and how models are specified,
etc. Available options are:

Sampling information mwg, global, block
Model definition moptimize, optimize, standalone
Model evaluator type d*, q*, e*, g*, v*
Other information fast

Thus, if the user is sampling from a previously formed model statement using optimize with a type d0 evaluator, and wished to sample in blocks the user
would specify amcmc_alginfo(A,("block","optimize","d0")).

amcmc_draw(A)

Use: Execute the sampler.

Setting up a problem as a structured object has advantages in
situations in which the user is interested in executing adaptive MCMC as a step in a larger algorithm, or in which the user must set up
a sequence of similar adaptive MCMC problems (or both). An additional set of commands which do not directly analogize with the use of the
function amcmc() follow:

amcmc_append(A, string scalar append)

Use: In the event that the user executes the command amcmc_draw(A) multiple times in sequence, by specifying append="append",
the results of the draw will be attached to the information on previous calls of amcmc_draw(A). Acceptance rates, the proposal distribution, etc.
are all updating accordingly. This is the default setting; unless the user specifies append="overwrite", all information about previous
draws will be retained. Specifying append="overwrite" does not restart the draw; everything is updated using past draws, but only the
most recent set of draws is retained. This is useful in large problems in which a sampler is used as a step in a larger sampling algorithm.

amcmc_reeval(A, string scalar reeval)

Use: In the event that the user executes the command amcmc_draw(A) as a step in larger algorithm,
the parameters that have been passed to the problem might change. It is then necessary for the algorithm to reevaluate the function
using the last drawn values with the new parameter values before proceeding. If the user wishes to do this, the user may simply
set reeval="reeval".

Step 3: Results

Once a run has been executed using amcmc_draw(), the user can access information about the draws, and also recover some information about
initial settings.

real matrix amcmc_results_draws(A)

Use: Returns the draws in a matrix form; each row represents a draw.

real colvector amcmc_results_vals(A)

Use: Returns the values of the proposal distribution corresponding with each draw.

real rowvector amcmc_results_arate(A)

Use: Returns acceptance rates; if a global scheme is used, this is a single value; otherwise, acceptance rates conform with the number of blocks
(for a block sampler) or the dimension of the target distribution (for a metropolis-within-gibbs sampler).

real scalar amcmc_results_passes(A)

Use: Returns the number of times the amcmc_draw() command has been issued.

real scalar amcmc_results_totaldraws(A)

Use: Returns the total number of draws.

real colvector amcmc_results_vals(A)

Use: Returns the values of the proposal distribution corresponding with each draw.

real matrix amcmc_results_acceptances(A)

Use: Returns a matrix of zeros and ones conformable with the sampling scheme, where a one marks whether or not a particular draw was
accepted, and a zero marks whether or not a draw was rejected. In short, the function returns the acceptance history of the draw.

real rowvector amcmc_results_propmean(A)

Use: Returns the mean of the proposal distribution at the end of the run.

real matrix amcmc_results_propvar(A)

Use: Returns the covariance matrix of the proposal distribution at the end of the run.

real rowvector amcmc_results_lastdraw(A)

Use: Returns only the last draw.

void amcmc_results_report(A)

Use: Produces a quick summary of a run, which includes the total passes,
draws per pass, the total number of draws, the average value of the target distribution,
the last value of the target distribution, the acceptance rate (average), and the burn-in period.

Additionally, one can obtain the values that the user specified by using amcmc_results_*(), where * can be any one of:
alginfo, noisy, blocks, damper, xinit, Vinit, or lambda.

Examples

Example 1: This example recasts example 2 above. To recap, the
interest lies in sampling from a multivariate normal mixture with dimension 2, where with probability p
the mean is m, and with probability p the mean is -m. p, m and Sigma - the covariance matrix of the distribution -
are passed as additional arguments of the function. A vector of zeros and an identity matrix are used as the starting values for the proposal
distribution:

real scalar mixnorm2(x,p,m1,m2,Sig)
>{
>dSig=1/sqrt(det(Sig))
>Siginv=invsym(Sig)
>val1=1/2*dSig*exp(-(x-m1)*Siginv*(x-m1)')
>val2=1/2*dSig*exp(-(x-m2)*Siginv*(x-m2)')
>return(ln(p*val1+(1-p)*val2))
>}

:p=1/2

:m1=1,1

:m2=-1,-1

:Sig=I(2):+.1

:Args=J(4,1,NULL)

:Args[1]=&p

:Args[2]=&m1

:Args[3]=&m2

:Args[4]=&Sig

:alginfo="standalone","global"

:A=amcmc_init()

:amcmc_alginfo(A,alginfo)

:amcmc_lnf(A,&mixnorm2())

:amcmc_draws(A,100000)

:amcmc_burn(A,10000)

:amcmc_args(A,Args)

:amcmc_xinit(A,J(1,2,0))

:amcmc_Vinit(A,I(2))

:amcmc_aopt(A,.34)

:amcmc_damper(A,1)

:amcmc_draw(A)

Example 2: Bayesian estimation of a mixed logit model. The ideas behind this example follow Chapter 12 of Train (2009).
The user is interested in estimating a logit model, but where parameters vary about a distribution at some group level.
There are three set of parameters: 1) beta_n, the values of the group parameters for each group, 2) b, the mean of the group
level parameters across groups, and 3) W, the covariance matrix of the group-level parameters.

Following the setup described in Train (2009), suppose that the prior distribution
for b is normal with arbitrarily large variance, while the prior on W is inverted
Wishart with K (the dimension of b) degrees of freedom and identity scale matrix. Under these conditions:

1. The distribution of b given beta_n and W is normal with mean=mean(beta_n), and covariance matrix W/N,
where N is the number of groups.

2. The distribution of W given b and beta_n is an inverted Wishart.

3. The distribution of beta_n given b and W is the product of a normal density
capturing the likelihood of the group's parameters given the mean and variance of group parameters across the sample,
and the likelihood of the individual's choices given beta_n.

The example begins by setting up the data, which is bangladesh.dta. The data set has information on
the use of contraceptives for individuals in different districts. The user posits that the coefficients vary
randomly by district in accordance with the model. Reading data into Mata:

. clear all

. webuse bangladesh
(NLS Women 14-24 in 1968)

. set seed 90210

. mata:

mata (type end to exit)

: st_view(X=.,.,"urban age child1 child2 child3")

: st_view(y=.,.,"c_use")

: st_view(id=.,.,"district")

: X=X,J(rows(X),1,1)

: m=panelsetup(id,1)

Mata now contains contraceptive use information in the vector y, explanatory variables (and a constant) in the vector X,
and an id code which is organized into a panel using panelsetup(). The next step is to code the functions
following Train's (2009) descriptions. For steps 1 and 2, a function producing draws from the
respective distributions is needed; step 3 will be set up to work with amcmc_*(), so
this function does not produce draws, but instead returns the log value of the parameter density conditional on data, choices, b, and W.
The notational convention in the necessary functions are that the mean and variance of group-level parameters are denoted by b and
W, the (matrix of) group-level parameters is denoted beta, while a set (really, a rowvector) of group-level parameters
are denoted by beta_n. The three required functions are:

: real matrix drawb_betaW(beta,W) {
>return(mean(beta)+rnormal(1,cols(beta),0,1)*cholesky(W)')
>}

: real matrix drawW_bbeta(beta,b)
>{
>v=rnormal(cols(b)+rows(beta),cols(b),0,1)
>S1=variance(beta:-b)
>S=invsym((cols(b)*I(cols(b))+rows(beta)*S1)/(cols(b)+rows(beta)))
>L=cholesky(S)
>R=(L*v')*(L*v')'/(cols(b)+rows(beta))
>return(invsym(R))
>}

: real scalar lnchoiceprob(beta_n,b,W,yn,Xn)
>{
>mus=rowsum(Xn:*beta_n)
>lnp=yn:*mus:-ln(1:+exp(mus))
>lnprior=-1/2*(beta_n-b)*invsym(W)*(beta_n-b)'-
> 1/2*ln(det(W))-cols(b)/2*ln(2*pi())
>return(sum(lnp)+lnprior)
>}

:

The next step is setting up a series of adaptive MCMC problems - one for each individual, using amcmc_*() and Mata's
J() function, which allows easy duplication of a single problem.

: Ap=amcmc_init()

: amcmc_damper(Ap,1)

: amcmc_arate(Ap,.4)

: amcmc_alginfo(Ap,("standalone","global"))

: amcmc_lnf(Ap,&lnchoiceprob())

: amcmc_draws(Ap,1)

: amcmc_append(Ap,"overwrite")

: amcmc_reeval(Ap,"reeval")

: A=J(rows(m),1,Ap)

:

Each problem is set up as a global drawing problem, where one draw is taken. In passing, the author's
experience is that it is sometimes helpful to
let each individual problem run a bit by specifying 5 or 10 draws, say, in this step for better mixing in the early stages of the
algorithm. The option "overwrite" specifies that
information on each draw is not to be stored as the algorithm proceeds. The option "reeval" specifies that
it is necessary to reevaluate the function, as the two arguments b and W are changed as the
algorithm proceeds. After some (pretty poor) starting values are
specified, an initial draw of the individual-level parameters is taken.
A loop can now be run to fill in the arguments. This is done by
setting up a pointer matrix, each row of which points to explanatory values and the mean and variance of the distribution. Now, initial
values for parameters are specified, and then all the separate amcmc problems are initiated in a loop.

: b=J(1,6,0)

: W=I(6)*6

:beta=b:+rnormal(rows(m),cols(X),0,1)*cholesky(W)'

: Args=J(rows(m),4,NULL)

: for (i=1;i<=rows(m);i++) {
>Args[i,1]=&b
>Args[i,2]=&W
>Args[i,3]=&panelsubmatrix(y,i,m)
>Args[i,4]=&panelsubmatrix(X,i,m)
>amcmc_args(A[i],Args[i,])
>amcmc_xinit(A[i],b)
>amcmc_Vinit(A[i],W)
>}

:

The algorithm is now implemented as follows (with 10000 total draws, and an initial value of individual
level parameters taken as a draw
from the normal distribution. The matrices bvals and Wvals are used to hold the draws of the mean and the covariance matrix:

: its=10000

: bvals=J(0,cols(beta),.)

: Wvals=J(0,cols(rowshape(W,1)),.)

: for (i=1;i<=its;i++) {
>b=drawb_betaW(beta,W/rows(m))
>W=drawW_bbeta(beta,b)
>bvals=bvals\b
>Wvals=Wvals\rowshape(W,1)
>for (j=1;j<=rows(m);j++) {
>amcmc_draw(A[j])
>beta[j,]=amcmc_results_lastdraw(A[j])
>}
>}

While in a typical application, one usually discards some number of initial values and thins results so as to
eliminate the autcorrelation inherent in MCMC sampling, here we just summarize the draws
for the mean of the parameter vector and the covariance matrix:

:mean(bvals)'
1
+
+
1 | .9653859143 |
2 | -.0363394328 |
3 | 1.399463261 |
4 | 1.696209143 |
5 | 1.656005711 |
6 | -2.121194287 |
+
+

:rowshape(mean(Wvals),6)
[symmetric]
1 2 3 4 5 6
+
+
1 | 1.692919191 |
2 | -.0236467397 .1201797953 |
3 | .7638987195 -.0098907509 1.356704652 |
4 | .5181705857 -.0080073073 .6239961521 1.122827391 |
5 | 1.087957638 -.0416998767 1.023219879 .6940770407 1.978052313 |
6 | -1.317031389 .0337024609 -1.044189251 -.812652926 -1.535208409 2.077123236 |
+
+

The example is in fact a fairly complete description as to how the user-written command bayesmlogit works.
For more information, type findit net bayesmlogit at the Stata command prompt. Further examples can be found in Baker (2013).

Additional Notes

amcmc() or amcmc_*() requires that the user install Ben Jann's help moremata set of commands.

Author

Matthew J. Baker, Hunter College and the Graduate Center, CUNY, matthew.baker@hunter.cuny.edu. Comments and suggestions are appreciated.

References

Andrieu, Christophe, and Johannes Thoms. 2008. A tutorial on adaptive MCMC. Statistics and Computing 18: 343-73.

Baker, Matthew J. 2013. Adaptive Markov chain Monte Carlo in Mata. Hunter College working paper 440.

Chernozukov, Victor, and Han Hong. 2003. An MCMC approach to classical estimation. Journal of Econometrics, 115(2): 293-346.

Train, Kenneth E. 2009. Discrete choice methods with simulation, 2nd. ed. Cambridge and New york: Cambridge University Press.

Also see

Online: help for mata