I am just a novice user, so there is surely a lot I do not know about the details of Stan. (I assume you meant to write sampling statement with the normal function and not N.) One difference between the approaches is that y ~ normal(mu, sigma) drops additive constants.
Quoting from the function reference:
y ~normal(mu, sigma)
Increment target log probability density with normal_lpdf(y | mu, sigma) dropping constant additive terms.
How exactly are the results different when you run your model?
Yeah @FJCC is correct. The difference is that the so-called “sampling notation” y ~ normal(mu, sigma) drops constants. The target += syntax keeps the constants and is also a lot more flexible, allowing you to use distributions for which we don’t have a built-in function.
Thanks for the reply Jonah, how about the inference parameters then, I see most of them are close to each other but not exactly the same. Is this something that can be explained by the normalization constant as well? When I ran the two inference, I used the same prior, same initialization point and same seed.
Good question. They should definitely give the same inferences for the parameters (within monte carlo error), but I also would have thought you’d get exactly the same values if, like you say, you also used the same inits and same seed. But I did a little experiment and I also found slight differences, so I think I’m also confused about something!
Maybe @bbbales2 or @rok_cesnovar knows what I’m missing here. In my example below, shouldn’t I get the exact same values of x or is that a faulty assumption?
# Fit two simple models, only differing on whether
# sampling notation or target+= is used
file1 <- write_stan_tempfile(
x ~ normal(1, 1);
file2 <- write_stan_tempfile(
target += normal_lpdf(x | 1, 1);
mod1 <- cmdstan_model(file1)
mod2 <- cmdstan_model(file2)
# initialize to 0 and use same seed
fit1 <- mod1$sample(init = 0, seed = 123)
fit2 <- mod2$sample(init = 0, seed = 123)
I expect different values for lp__ but I also get slightly different results for x, certainly within monte carlo error, but it still surprised me:
> fit1$print("x", "mean", "sd")
variable mean sd
x 1.01 0.99
> fit2$print("x", "mean", "sd")
variable mean sd
x 1.00 1.01
I notice that the number of iterations in @LeonCaesa’s example is only 50, Rhat is almost always far above the recommended 1.01 (is that the latest recommendation?) upper limit, and n_eff is small. I would try running with iter = 2000, if that does not take too long.
Yeah, this is probably it. The large Rhat indicates that Stan isn’t mixing and the results are not to be relied upon (so it’s not too surprising you saw differences). You should run it until the Rhats are small before you do comparisons.