In the simple case of normally distributed data with unknown mean and variance, Jeffrey's prior is given by $$p(\mu, \sigma^2)=\frac{1}{\sigma^2}.$$ How can I define such a prior in the Stan language, i.e. how do I need to change the model statement below to obtain the desired result? (The model statement now is for the prior $p(\mu, \sigma^2)\propto 1.$

data {
 int<lower=0> n; // obs in group x
 real x[n];
parameters {
 real muX;
 real<lower=0> sigmaSquared;
 real postPred;
transformed parameters
 real<lower=0> sigmaX;
 sigmaX <- sqrt(sigmaSquared);
model {
 x ~ normal(muX, sigmaX);
 postPred ~ normal(muX, sigmaX);
  • 1
    $\begingroup$ You can't (at least I do not know any way to) directly define improper prior other than the default priors (the one you mentioned). You could instead use zero-ones trick to define an entire likelihood or, in your case, likelihood + prior. This should be very similar to the case of OpenBUGS software. $\endgroup$ – Tomas Jun 12 '15 at 19:24
  • $\begingroup$ Thanks. Could you elaborate on what you mean by 'zero-ones trick' and where I can read up on how to define an entire likelihood? $\endgroup$ – Nussig Jun 12 '15 at 19:26
  • $\begingroup$ In fact you don't need zero-one trick for the entire likelihood, only for your new prior. You can find info on how to construct them in OpenBUGS here: users.aims.ac.za/~mackay/BUGS/Manuals/Tricks.html Under the "Specifying a new prior distribution". Theta would be your $\sigma$ in that specific case $\endgroup$ – Tomas Jun 12 '15 at 19:33

You can define a proper or improper prior in the Stan language using the increment_log_prob() function, which will add its input to the accumulated log-posterior value that is used in the Metropolis step to decide whether to accept or reject a proposal for the parameters.

In your example, the model block would need to include the new line increment_log_block(-log(sigmaSquared)); However, some people (e.g. Jaynes) argue that the Jeffreys prior is only appropriate for scale parameters, in which case you could reparameterize your model in terms of the standard deviation (sigmaX) rather than the variance (sigmaSquared). Also, what I assume is your attempt to draw from the posterior predictive distribution of x should be in a generated quantities block. Putting all three pieces together, it would look like: data { int<lower=0> n; // obs in group x real x[n]; } parameters { real muX; real<lower=0> sigmaX; } model { x ~ normal(muX, sigmaX); increment_log_prob(-log(sigmaX)); } generated quantities { real postPred; postPred <- normal_rng(muX, sigmaX); }

  • $\begingroup$ Thank you, Ben Goodrich, for your comment and your edit. Could you please elaborate (or provide a source) on the difference of having postPred in the generated quantities instead of in the model statement? $\endgroup$ – Nussig Jun 13 '15 at 17:23
  • $\begingroup$ Anything that can go into generated quantities should go into generated quantities block, which is evaluated with plain double precision values. Anything in parameters, transformed parameters, or model involves Stan's custom scalar types, which calculate the gradient via autodifferentiation and hence all operations involving them are slower. Conceptually, I think the code is easier to understand if you draw from the posterior distribution of the parameters and the posterior predictive distribution of the generated quantities. $\endgroup$ – Ben Goodrich Jun 13 '15 at 18:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.