Estimating step_size in diag_nuts

[Yuanjun Gao]

Hi. I have two quick questions about estimating step_sizes in diag_nuts.hpp

1. When estimating the step_sizes, why do you do the standardization(by using "normalizer") to make the square sum of the step_sizes equal to the number of parameters? The relevant code is

if (step_size_sq_sum > 0.0) {
_x_sum_n = 0;
double normalizer = sqrt((double)_step_sizes.size())
/ sqrt(step_size_sq_sum);
for (size_t i = 0; i < _step_sizes.size(); i++)
_step_sizes[i] *= normalizer;
}

2. In each iteration, when calculating _x_sum and _xsq_sum. You are not only using the final result of the parameters. Instead, you use "nvalid" to judge whether or not to add the "newx", (as the annotation says, you "update running statistics if point is in slice") such that during each iteration, you may get several "newx". Why do you do this? The relevant code is

if (nvalid) {
_x_sum_n++;
for (size_t i = 0; i < newx.size(); i++) {
_x_sum[i] += newx[i];
_xsq_sum[i] += newx[i] * newx[i];
}
}

Thank you for your help!

[Bob Carpenter]

On 12/3/12 8:48 AM, Yuanjun Gao wrote:
> Hi. I have two quick questions about estimating step_sizes in diag_nuts.hpp

Good questions! This is helpful for me to think about, too.

> 1. When estimating the step_sizes, why do you do the standardization(by using "normalizer") to make the square sum of the step_sizes equal to the number of parameters? The relevant code is
>
> if (step_size_sq_sum > 0.0) {
> _x_sum_n = 0;
> double normalizer = sqrt((double)_step_sizes.size())
> / sqrt(step_size_sq_sum);
> for (size_t i = 0; i < _step_sizes.size(); i++)
> _step_sizes[i] *= normalizer;
> }
>

I'm pretty sure this is because the overall step size
is also being adapted --- think of it as an identifiability
issue.

> 2. In each iteration, when calculating _x_sum and _xsq_sum. You are not only using the final result of the parameters. Instead, you use "nvalid" to judge whether or not to add the "newx", (as the annotation says, you "update running statistics if point is in slice") such that during each iteration, you may get several "newx". Why do you do this? The relevant code is
>
> if (nvalid) {
> _x_sum_n++;
> for (size_t i = 0; i < newx.size(); i++) {
> _x_sum[i] += newx[i];
> _xsq_sum[i] += newx[i] * newx[i];
> }
> }

That will use adaptation all along the path of the particle,
which should be more effective an estimator than just using the
endpoint. If you tried to include points that are out of the slice,
you could wind up with not-a-number (NaN) values. I'm also
guessing that including points in the support but with low
probabilities might also reduce the efficiency of the estimator.

I'm curious what Matt has to say.

- Bob

[Matt Hoffman]

Oops, I just replied to Yuanjun instead of to the list. Here's my
response, which more or less echoes Bob's:

1. There is a sort of an identifiability problem with the global step
size epsilon and the vector of per-variable step sizes—one could
multiply the global step size by x and divide all of the local step
sizes by x and get the same result. If the overall magnitude of the
step_sizes vector changed dramatically, this would throw off the
estimation of the global step size parameter. So it's more stable if
we ensure that only the relative weights in the step size vector
change, and let the overall step size be controlled by epsilon.

2. There is useful information throughout each trajectory. Any
in-slice point is (more or less) a valid sample from the posterior, so
there's no reason to throw away that information when estimating the
relative variances of each dimension. So that's why I use all of the
in-slice points on the trajectory, rather than just the point that's
ultimately sampled. This removes one source of noise when estimating
the variances.

Hope that helps!

Matt

> --
>
>