diffuse inits -> low prob -> high gradients -> NaN

[Bob Carpenter]

(Marcus: yes, you're absolutely right -- it's a model
problem. And we'll never solve all of those. You can
specify inits now in the command, but there's no good way
yet to specify them programatically -- it's on the to-do
list and we may set the summer student to writing and
testing the random generators we'll need.)

(Jiqiang: you still only need a single seed;
we took the approach Marcus suggested.)

So here's the story:

With init retries, the BUGS vol1/salm model at least took off.
But it crashed after 12 iterations with NaN everywhere.

What appears to be happening (see below for a dump)
is that random inits in (-2,2) on the unconstrained scale can
lead to very small initial log probabilities (-1E34) and
very large gradients (1E34) even when the initial log prob
isn't NaN. Finite diffs fail by returning zero (see below),
but the actual model gradients are all finite.

In contrast, all-zero inits work fine, with all-zero inits,
the log prob is -2E1 (yes, E1), and the biggest gradient is 1E5,
the maximum error from finite diffs is 1E-8 and everything runs
fine.

Daniel also brought up the problem that the random momenta
on the unit normal scale can cause problems, but those should
only lead to rejections if they're really badly tuned. And
the varying step size (equivalently variable diagonal mass matrix)
should help.

I think we either need to

a) more agressively lower the step size in the sampler
during warmup (but then it has to be allowed to creep back
up when we get into higher probability regions), or

b) trim the init interval from (-2,2) to something much
less agressive, say (-.25,.25).

We'll still run into problems if all-zero inits don't behave
well because we're initializing around them. For the salm model,
the predictor vector x has values ranging from 0 to 1000.
The salm.stan model centers them, but doesn't rescale.
Presumably the Poisson's blowing out in exponentiating
the linear basis in the model (y[i,j] ~ poisson(exp(...)));
I don't see how to fix that by taking a log(lambda) as the
parameter instead of lambda -- you wind up needing to
exponentiate to calc the log probability anyway.

We'll have an issue if someone initializes a precision value
uniformly, as in

tau ~ uniform(0,1000);
sigma <- 1 / sqrt(tau);

In this case, a 0-init on tau leads to 1/sqrt(500) = .04.
On the other hand,

sigma ~ uniform(0,1000);

will start with sigma = 500, which is probably better
behaved than having the scale too low.

Of course, there's no way to solve our init problems
in general. An adversarial model designer will still
be able to trash any general init scheme.

Here's the example from the salm model, with some ill-behaved
random inits:

-----------------------------
TEST GRADIENT MODE

Log probability=-1.23088e+34

param idx value model finite diff error
0 -0.779718 -1.23088e+34 -1.23088e+34 3.02684e+25
1 -1.01175 2.56122e+34 2.56122e+34 3.08556e+25
2 -0.306507 3.0216e+36 3.0216e+36 -3.02353e+28
3 0.477941 -12.5309 0 -12.5309
4 -0.221135 -1.74369e+33 -1.74369e+33 -4.42628e+24
5 1.37261 -8.5827e+33 -8.5827e+33 -2.20417e+25
6 -0.480064 -1.34592e+33 -1.34592e+33 -3.67377e+24
7 1.77561 -2.97122e+32 -2.97122e+32 -3.64414e+23
8 1.70499 -2.76864e+32 -2.76864e+32 1.11293e+24
9 0.213972 -6.2334e+31 -6.2334e+31 2.89479e+23
10 0.209229 -2.48138e+28 -2.48143e+28 5.01093e+23
11 0.731258 -4.18224e+28 -4.18222e+28 -1.23845e+23
12 1.70378 -1.10604e+29 -1.10603e+29 -3.88902e+23
13 -1.61234 -1.87167e+18 0 -1.87167e+18
14 1.53896 -4.37342e+19 0 -4.37342e+19
15 -0.467814 -5.87883e+18 0 -5.87883e+18
16 -0.467691 33.7543 0 33.7543
17 -1.18276 39.9075 0 39.9075
18 -1.96958 44.1764 0 44.1764
19 1.07951 18.259 0 18.259
20 -0.675662 28.0897 0 28.0897
21 -1.62609 44.6225 0 44.6225
-------------------------------

(The 0 values for finite diff were initially worrying, but they match
just fine when the log probability is in the neighborhood of -1E7
rather than -1E34.)

- Bob

[Matt Hoffman]

Hi all,

This was due to a bug in nuts_diag.hpp where it would estimate 0
variance for every dimension, leading to per-variable step sizes of
0/0=NaN. I'll push the fix later today.

Matt