Re: Something wrong in the code?

[Bob Carpenter]

[discussion moved to stan-dev list -- if you
don't have posting priveleges yet, Yuanjun, we
can give them to you]

On 11/14/12 11:45 PM, Yuanjun Gao wrote:
> Dear Bob and Daniel,
>
> While I am reading the implementation of nuts_diag, I found something weird for me and I think there might be some flaws
> in the implementation. So I hope you can help me to check.
>
> First, in stan/mcmc/util.hpp, the implementation of "rescaled_leapfrog" seems strange. based on R.Neal's paper
> (2.28)-(2.30), the step_sizes(which is just the inverse mass), should only appear when calculating the new position, but
> shouldn't appear when calculating the new momentum. So I think for the 89th and 99th line, the step_sizes[i] shouldn't
> appear. Right? I paste the code in until.hpp as follows
>
> 86 for (size_t i = 0; i < m.size(); i++)
> 87 m[i] += 0.5 * epsilon * step_sizes[i] * g[i];
> 88 for (size_t i = 0; i < x.size(); i++)
> 89 x[i] += epsilon * step_sizes[i] * m[i];
> 90 double logp;
> 91 try {
> 92 logp = model.grad_log_prob
> <http://mc-stan.org/api/classstan_1_1model_1_1prob__grad.html#ab460008516b71b39c546526f5faf710b>(x, z, g, output_msgs);
> 93 } catch (std::domain_error e <http://mc-stan.org/api/namespacestan_1_1math.html#a3b650a0131d41167ef4837ecc7d02be5>) {
> 94 write_error_msgs <http://mc-stan.org/api/namespacestan_1_1mcmc.html#a50d1d353df9107598d29b6d79f10cd91>(error_msgs,e);
> 95 logp = -std::numeric_limits<double>::infinity();
> 96 }
> 97 for (size_t i = 0; i < m.size(); i++)
> 98 m[i] += 0.5 * epsilon * step_sizes[i] * g[i];
> 99 return logp;

It's a varying step size implementation, not a varying
mass one. See section 5.4.2.4 of Radford's
paper for a discussion.

Also, back of the envelope, we have a lot of trust
in the way it is now because (1) we've tried it on
a lot of models and verified answers against BUGS,
(2) if it were wrong, we'd see lots of Metropolis
rejections because the Hamiltonian would diverge, and
(3) Matt wrote it.

> Secondly, in "nuts_diag.hpp", in the "build_tree" function, the program calculate newH (the new hamilton quantity) as:
>
> 360 double newH = newlogp - 0.5 * stan::math::dot_self
> <http://mc-stan.org/api/namespacestan_1_1math.html#ac81ab85d390fcb594884773a85fac657>(minus);
>
> But now since the mass matrix is not an identity but a diagonal matrix, we should consider the mass as well, so in this
> case, we should use:
>
> newH = newlogp - 0.5 * sum( mminus[i]^2 / mass[i] ) ; (sorry for being too lazy to translate it into C++�)
>
> I am looking forward to your reply. Thank you very much!

Again, we're just using varying step sizes.

My understanding of the plan for estimating the mass matrix
is to have a single step size and a mass matrix. If we do that,
then you need to change the kinetic energy term in the Hamiltonian
and the updates as you suggest. For the kinetc, it's now

kinetic energy = p' * inv(M) * p

I also believe the plan for warmup is as follows:

draw 1--20 with simple step size adaptation
use draws 1--20 to estimate M

draw 21--40 with step size adaptation and mass matrix M
use draws 21--40 for new estimate of M
(i.e., discard draws 1--20 for the estimate of M)

draw 41--80 and use these draws for new estimate of M

....

Then after warmup, mass matrix M and step size remain fixed.

What I was suggesting as an implementation strategy was to
reimplement a new sampler based on a mass matrix instead of
varying step sizes, then recast what we're doing now as a diagonal
mass matrix, then finally go to a full-blown mass matrix.

- Bob

[Bob Carpenter]

On 11/15/12 12:58 AM, Yuanjun Gao wrote:
> Dear Bob,
>
> Thank you for your reply, but which paper of Radford are your referring to? I looked at <MCMC using Hamiltonian dynamics>(2010) but just found section 5.4, where is 5.4.2.4?

The MCMC Handbook paper:

http://www.mcmchandbook.net/HandbookSampleChapters.html

> BTW, I don't know what is stan-dev list�

It's the mailing list for Stan developers. So we
all get it.

Here's the Stan Google Code page for Stan:

https://code.google.com/p/stan/

And here's the link to the stan-dev list:

https://groups.google.com/forum/?fromgroups#!forum/stan-dev

You need to request to join and we need to give you permission
to post.

- Bob