centered vs non-centered parameters: very different optimization results expected?

[Charles Driver]

I was surprised that in the model I've been working on, whenever I optimized it, the hyper variances ran out to the upper bound very quickly. To understand what was going on I set up a toy in R, which indeed shows that while point estimates of centered hyper variances are in general too low, those of noncentered variances appear unbounded in their desire to be large...

Is there some reason I should have expected this? Have I misunderstood something crucial about the noncentered approach, or about optimizing in Stan?  Below is the toy script, thanks!


hn=50
hmu<-2
hsd<-.2
ar=.6
mu<-rnorm(hn,hmu,hsd)

n=10
diffusion=.4
x<-matrix(rnorm(hn,1/(1-ar)*mu,sqrt(diffusion^2/(1-ar^2))),ncol=hn)
for(i in 2:n){
  x <- rbind(x,x[i-1,]*ar+mu+rnorm(hn,0,diffusion))
}


require(rstan)
noncentered<-"
data{
int hn;
int n;
vector[n] x[hn];
}
parameters{
real<lower=0.0001,upper=1> ar;
real<lower=.001,upper=10> diffusion;
real<lower=.001,upper=10> hsd;
real hmu;
vector[hn] mudev;
}
model{
vector[hn] mu;
mudev~normal(0,1);
mu<-hmu+mudev*hsd;
for(subjecti in 1:hn){
x[subjecti,1]~normal(1/(1-ar)*mu[subjecti],((diffusion^2)/(1-ar^2)));
for(rowi in 2:n) x[subjecti,rowi]~normal(mu[subjecti]+ar*x[subjecti,rowi-1],diffusion^2);
}
}"



centered<-"
data{
int hn;
int n;
vector[n] x[hn];
}
parameters{
real<lower=0.0001,upper=1> ar;
real<lower=.001,upper=10> diffusion;
real<lower=.001,upper=10> hsd;
real hmu;
vector[hn] mu;
}
model{
mu~normal(hmu,hsd);
for(subjecti in 1:hn){
x[subjecti,1]~normal(1/(1-ar)*mu[subjecti],((diffusion^2)/(1-ar^2)));
for(rowi in 2:n) x[subjecti,rowi]~normal(mu[subjecti]+ar*x[subjecti,rowi-1],diffusion^2);
}
}"


standata<-list(x=t(x),hn=hn,n=n)

centerednuts<-stan(model_code = c(centered),
  data = standata, iter = 1000, chains = 3,
  cores=3)
centerednuts

noncenterednuts<-stan(model_code = c(noncentered),
  data = standata, iter = 1000, chains = 3,
  cores=3)
noncenterednuts

centeredoptim<- optimizing(object = centerednuts@stanmodel,
  algorithm='BFGS',
  as_vector=F,
  tol_obj=1e-10, tol_grad=1e-10,tol_param=1e-10,tol_rel_grad=0, tol_rel_obj=0,
  data = standata)
centeredoptim

noncenteredoptim<- optimizing(object = noncenterednuts@stanmodel,
  algorithm='BFGS',
  as_vector=F,
  tol_obj=1e-10, tol_grad=1e-10,tol_param=1e-10,tol_rel_grad=0, tol_rel_obj=0,
  data = standata)
noncenteredoptim

[Bob Carpenter]

That should give you the same answer as

> vector[hn] mu;
> mudev~normal(0,1);
> mu<-hmu+mudev*hsd;

mu ~ normal(hmu, hsd);

assuming you declare mu as an unconstrained parameter.

You shouldn't need that lower=0.001 for scales in Stan. If you
do need it, it's indicating a problem with your model and will
lead to computational problems trying to sample near the boundary.
We also don't recommend upper-bounding scale parameters.

- Bob

> On May 10, 2016, at 2:14 PM, Charles Driver <cdriv...@gmail.com> wrote:
>
> I was surprised that in the model I've been working on, whenever I optimized it, the hyper variances ran out to the upper bound very quickly. To understand what was going on I set up a toy in R, which indeed shows that while point estimates of centered hyper variances are in general too low, those of noncentered variances appear unbounded in their desire to be large...
>
> Is there some reason I should have expected this? Have I misunderstood something crucial about the noncentered approach, or about optimizing in Stan? Below is the toy script, thanks!
>
>
> hn=50
> hmu<-2
> hsd<-.2
> ar=.6
> mu<-rnorm(hn,hmu,hsd)
>
> n=10
> diffusion=.4

> x<-matrix(rnorm(hn,1/(1-ar)*hx,sqrt(diffusion^2/(1-ar^2))),ncol=hn)
> for(i in 2:n){
> x <- rbind(x,x[i-1,]*ar+hx+rnorm(hn,0,diffusion))

> --
> 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.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Charles Driver]

Right, I also expected that should give the same result, it just, doesn't :) Mu is unconstrained. I went back and unconstrained hsd for the sake of it, now in the noncentered approach it reaches crazy high or low values while optimizing.

[Michael Betancourt]

It gives the same posterior expectations but will

yield different modes as they are densities with

respect to different parameters.  Of course the

mode is a terrible and shouldn’t be used, in

which case this sensitivity is irrelevant.

On May 12, 2016, at 5:06 PM, Charles Driver <cdriv...@gmail.com> wrote:

Right, I also expected that should give the same result, it just, doesn't :) Mu is unconstrained. I went back and unconstrained hmu for the sake of it, now in the noncentered approach it reaches crazy high or low values while optimizing.

[Bob Carpenter]

There's no Jacobian involved here, so why would they
have different modes (mu, not mudev, that is)?

- Bob

[Michael Betancourt]

Because you’re optimizing with respect to theta_tilde, not
theta. It’s a different parameter space (i.e. the centered
and non-centered parameterizations aren’t different
parameterizations of the same system, they are
fundamentally different spaces that can be manipulated
into the same system).

[Bob Carpenter]

Surprising (the result, not that Michael's right about calc!).
It looks like we're optimizing f in one case and h in the second,
where

f(y_raw, mu, sigma,theta) = g(sigma * y_raw + mu, mu, sigma, theta)
+ log normal(y_raw | 0, 1)

vs.

h(y, mu, sigma, theta) = g(y, mu, sigma, theta)
+ log normal(y | mu, sigma)

= g(y, mu, sigma, theta)
+ log normal((y - mu)/sigma | 0, 1)
- log(sigma)

where g is just the rest of the model log density.
The only difference is the -log(sigma) term, but that will
cause different densities.

If I did the algebra right, this should go in the manual.

- Bob

[Krzysztof Sakrejda]

+1 for going in the manual! Krzysztof

>>>> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

>>>> To post to this group, send email to stan-...@googlegroups.com.
>>>> For more options, visit https://groups.google.com/d/optout.
>>>
>>>
>>> --
>>> 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+unsubscribe@googlegroups.com.

>>> To post to this group, send email to stan-...@googlegroups.com.
>>> For more options, visit https://groups.google.com/d/optout.
>>
>> --
>> 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+unsubscribe@googlegroups.com.

>> To post to this group, send email to stan-...@googlegroups.com.
>> For more options, visit https://groups.google.com/d/optout.
>
> --
> 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+unsubscribe@googlegroups.com.

[Charles Driver]

Thanks for the input all. It's still not totally clear to me how I can increment the log prob to account for different parameterizations, but it's not so important - just good to understand and confirm that it's expected!

[Bob Carpenter]

> On May 19, 2016, at 12:22 PM, Charles Driver <cdriv...@gmail.com> wrote:
>
> Thanks for the input all. It's still not totally clear to me how I can increment the log prob to account for different parameterizations, but it's not so important

I have to disagree --- it's critical to understand what
you're computing and thus to make sure you're computing
what you think you're computing.

I wrote up a case study of a very simple example
here:

http://mc-stan.org/documentation/case-studies/mle-params.html

- Bob

[Charles Driver]

The write up is nice, thanks. And I agree it's important to 'know what your computing', but there are different levels of that understanding, sometimes a conceptual understanding is an adequate state, given a mix of competing worldly constraints!