Differences between ~normal() and target += normal_lpdf()


#1

Hi - I am writing up some examples, and am encountering some strange behavior. I have the following model straight from the manual:

data{
  int N;
  int K;
  matrix[N, K] X;
  vector[N] y;
}

parameters{
  real alpha;
  vector[K] beta;
  real<lower = 0> sigma;
}

model{
  y ~ normal(X * beta + alpha, sigma);
}

This give the following estimates:

          mean se_mean   sd     5%    95% n_eff Rhat
alpha    40.92    0.10 2.84  36.38  45.59   841    1
beta[1]  -1.30    0.02 0.68  -2.40  -0.17   913    1
beta[2]   0.01    0.00 0.01  -0.01   0.03   791    1
beta[3]  -0.02    0.00 0.01  -0.04   0.00  1115    1
beta[4]  -3.88    0.03 1.05  -5.60  -2.08   924    1
sigma     2.62    0.01 0.38   2.08   3.31   847    1
lp__    -45.22    0.08 1.95 -48.83 -42.81   539    1

Changing the likelihood line to:

target += normal_lpdf(y | X * beta + alpha, sigma);

Results in these estimates:

           mean se_mean   sd     5%    95% n_eff Rhat
alpha     33.78    0.11 2.97  28.55  38.40   668 1.01
betas[1]  -0.29    0.02 0.71  -1.42   0.90   980 1.00
betas[2]  -0.01    0.00 0.01  -0.04   0.01   764 1.00
betas[3]  -0.02    0.00 0.01  -0.04   0.00  1178 1.00
betas[4]  -1.87    0.04 1.05  -3.48  -0.05   834 1.00
sigma      2.91    0.02 0.48   2.24   3.81   709 1.00
lp__     -91.35    0.09 2.03 -95.39 -88.80   562 1.00

The data are the mtcars dataset in R, with the model mpg ~ cyl + disp + hp + wt.

This happens in both rstan and cmdstan. Any help appreciated.


#2

They both result in the same density wrt to the sampling. The difference is that normal_lpdf returns the correctly normalized density. It is thus a tiny bit slower as the normalizing constants are being calculated each time. The ~ normal notation does drop the normalizing constants.

So within MC error your results should match (and they do as it looks like).


#3

From the example above, if the mc standard error is under se_mean, it looks like the alpha's are quite significantly different, with a difference of 7 and a mc error of 0.1?


#4

OK - I thought I was misunderstanding the differences between these two, but given your response this might be a bug.

If you look closely, the two samples do not match. The mean for beta[4] in the first one is outside of the interval for the second. The mean for beta[1] on the first is very close to the edge of the interval for the second. With an effective size ~800, there should not be this large of a difference.

For reference, the ~normal() notation gives the same results as rstanarm::stan_glm().

My setup:

Linux 4.15.7
gcc 7.3.0
cmdstan 2.17.1
R 3.4.3
rstan 2.17.3


#5

@mespe Can I convince you to put a minimal reproducible example in the relevant repo (whatever interface you’re running)?


#6

I would be happy to, though I observed it in multiple interfaces (cmdstan and rstan). Should it go to https://github.com/stan-dev/stan?


#7

Sort of doesn’t matter because it’ll end up getting moved to either stan-dev/stan or stan-dev/math

Somebody just needs to run multiple chains, confirm that data is passed in correctly, confirm that the two are different, confirm that it’s not a mixing problem from the model, check what c++ is generated, etc…

if there are actually problems with the normal log density we should hand out prizes for finding them :)


#8

Issue opened on github for stan. Just let me know if you need additional pieces of information.


#9

Oh… you are right. I did not look that closely. Have you

  • Run the above with different seeds?
  • Rerun the thing with at least weakly informative priors?
  • Can you check the sampler diagnostics? (stepsize, masses)… if these terms don’t mean too much to you, don’t mind.

I would suspect that there is something odd going on wrt to the warmup and no priors are not a good thing - but you are right in that results should not differ that much.


#10

I have tried with and without priors of varying strength, different starting seeds, different number of iterations. I pulled the priors off the example model to make it the bare minimum model necessary to see the issue.

On my machine, it is pretty consistent. That is why I thought I was not understanding the difference between these two fully.

I can put together a more comprehensive test case if that is helpful.


#11

I just ran this on a different machine and did not experience any differences, but still get those differences on my laptop.

I need to dig into what happened on this laptop to cause this. Any ideas where to start?


#12

that’s good news!

You should wipe out all binaries, then check your compiler settings and rerun on your laptop.


#13

After few hours, I think I figured it out. My kernel was updated, but my laptop was not restarted.

A reboot later, both are now the same to MCMC error. Sigh…


#14

This seems to be a common problem people are running into. Seems like the OS should suggest a restart after updating!