Initialization not working as expected

[Brian Lau]

I came across some confusing behavior in initialization (CmdStan 2.1.0, OSX 10.7.5, g++-4.2).

Using the example bernoulli model from the manual (in src/models/basic_estimators)

./bernoulli sample save_warmup=1 data file=bernoulli.data.R

Produces a random initial value for theta as expected by the default init=2

./bernoulli sample save_warmup=1 init=0 data file=bernoulli.data.R

Results in the initial value for theta = 0.5. Shouldn't this be zero?

./bernoulli sample save_warmup=1 init=init.R data file=bernoulli.data.R

where init.R is a dump file with theta <- 0 produces the message "Rejecting user-specified initialization because of vanishing density." I guess this happens when an init exactly equals a bound? This call also returns an exit value of 0. Maybe it should return a more informative exit value, like when a initial value exceeds a bound. Changing init.R to theta <- 0.75 works and produces an expected initial value of 0.75.

It got even stranger when I tested a model with a vector parameter. I tried a very simple regression (files attached) with the following:

./mlr sample save_warmup=1 init=0 data file=data.R

which fails to produce initial values equal to zero.

./mlr sample save_warmup=1 init=init.R data file=data.R

also usually does not produce the correct initial values. Close but not quite. I say usually because I've found that when sampling from more than one chain, I occasionally get a chain with exactly the correct initial conditions, for all parameters.

I'm assuming the first row of samples in the csv files represent the initial conditions. Is this correct?

Thanks for any insights,

-b

[Bob Carpenter]

On 1/18/14, 2:15 PM, Brian Lau wrote:
> I came across some confusing behavior in initialization (CmdStan 2.1.0, OSX 10.7.5, g++-4.2).
>
> Using the example bernoulli model from the manual (in src/models/basic_estimators)
>
> ./bernoulli sample save_warmup=1 data file=bernoulli.data.R
>
> Produces a random initial value for theta as expected by the default init=2
>
> ./bernoulli sample save_warmup=1 init=0 data file=bernoulli.data.R
>
> Results in the initial value for theta = 0.5. Shouldn't this be zero?

0.5 is correct. It's 0 on the *unconstrained* scale, which is then
inverse-logit transformed to (0,1) resulting in 0.5.
The [-2,2] uniform init is also on the unconstrained scale.
The reason we do this is that we're then guaranteed to get
a valid initialization point if the model's parameters are declared
with appropriate constraints (so support is over all valid variables).

The init files, on the other hand, are on the constrained scale,
and it's up to the user to make sure they satisfy the constraints.

> ./bernoulli sample save_warmup=1 init=init.R data file=bernoulli.data.R
>
> where init.R is a dump file with theta <- 0 produces the message "Rejecting user-specified initialization because of
> vanishing density." I guess this happens when an init exactly equals a bound? This call also returns an exit value of 0.
> Maybe it should return a more informative exit value, like when a initial value exceeds a bound. Changing init.R to
> theta <- 0.75 works and produces an expected initial value of 0.75.
>
> It got even stranger when I tested a model with a vector parameter. I tried a very simple regression (files attached)
> with the following:
>
> ./mlr sample save_warmup=1 init=0 data file=data.R
>
> which fails to produce initial values equal to zero.
>
> ./mlr sample save_warmup=1 init=init.R data file=data.R
>
> also usually does not produce the correct initial values. Close but not quite. I say usually because I've found that
> when sampling from more than one chain, I occasionally get a chain with exactly the correct initial conditions, for all
> parameters.

> I'm assuming the first row of samples in the csv files represent the initial conditions. Is this correct?

Yes, that's supposed to happen if save_warmup=1.

You've definitely tickled a bug here. If I do this:

~/temp2/lau$ ./mlr sample save_warmup=1 init=init.R data file=data.R

it saves the initialization as the first sample, whereas if I do this:

./mlr sample num_samples=100 num_warmup=100 save_warmup=1 init=init.R data file=data.R

it saves 200 iteration values, but the first one's not the initialization.

- Bob

[Andrew Gelman]

We've changed our default to [-.5,.5], right?

[Bob Carpenter]

Not yet. More on stan-dev in a second.

- Bob

[Bob Carpenter]

I just checked with Michael and so far, there's no guarantee
that the initial values will be saved as the first iteration
if the num_warmup or num_samples is set.

We do plan to fix it so that the init is always the first
iteration, but it's not going to happen for the upcoming
2.1.1 release, which will fix a number of outstanding issues.
Maybe in 2.2.0.

- Bob

[Steve]

Just to clarify - the initial values in the init.R file are in fact used as the initial values on the constrained scale, and the first values saved are the 2nd set of values in the chain?

[Bob Carpenter]

On 1/19/14, 3:02 AM, Steve wrote:
> Just to clarify - the initial values in the init.R file are in fact used as the initial values on the constrained scale,
> and the first values saved are the 2nd set of values in the chain?

Right, but the "2nd set of values" notion is a bit subtle because
there's an adaptation step at the get go that looks for a suitable
step size to make progress.

The first phase of adaptation (after finding a workable step size)
then uses that initial step size with a unit mass matrix
(same step size in every dimension) in an attempt to find the high
mass volume of the posterior.

The second phase of adaptation then alternatively looks for a step
size that leads to the right rejection rate (set at 0.8 in Stan 2.1.0
byd efault) and estimates the mass matrix (effectively
the scale of each parameter in the default setting, which is a diagonal
mass matrix).

The third and final phase of adaptation then fine tunes the step size
with the final estimate of the mass matrix from the second adaptation phase.

The length of all of these phases is tunable in Stan 2.1.0, but it's
probably easier to just run for enough iterations that the chains
start mixing well after adaptation.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.

[Brian Lau]

Thanks for all the info Bob.

 

You've definitely tickled a bug here.  If I do this:

   ~/temp2/lau$ ./mlr sample save_warmup=1 init=init.R data file=data.R

it saves the initialization as the first sample, whereas if I do this:

   ./mlr sample num_samples=100 num_warmup=100 save_warmup=1 init=init.R data file=data.R

it saves 200 iteration values, but the first one's not the initialization.

FYI, neither of these commands give me the initialization as the first sample. Example output from the first command above attached in case its useful.

-b

[Bob Carpenter]

We took the discussion up on the dev list.

The current state of Stan is that the first line output
if save_warmup=1 is the first iteration *after* the
inits. It may be close to the inits if the sampler didn't
move far or exactly the inits if the first move was rejected.

We're going to fix this in a future release so that we
always save the init value as the first iteration. This will
mainly be useful for random inits, because there's no way to
recover the initialized values otherwise.

- Bob

[Steve]

Thanks Bob. All very helpful. The proposed change will also be helpful as a verification to those who supply initial values that their syntax was correct and STAN actually used the given init values.

[Bob Carpenter]

On 1/19/14, 6:10 PM, Steve wrote:
> Thanks Bob. All very helpful. The proposed change

[always including initial values in output]

> will also be helpful as a verification to those who supply initial
> values that their syntax was correct and STAN actually used the given init values.

Good point on verification --- the more ways we have for
users (including us!) to verify the model's working, the
better.

For now you can also print the values of parameters
in the model. It will get called first for the initial values.

- Bob