I'm just learning Bayesian statistics with stan, so bear with me.

Let's say we have one observed data points:

tweets_dat <- list(J = 1,
                    y = 2500,
                    sigma = ?)

I've left the error term incomplete as I don't yet know exactly what should go there.

I know that the mean from my observations should have mean of 1100. In other words, if the 2500 is quite high. I estimate it is two standard deviations away from the mean.

How can I make my model reflect this prior knowledge?

I'm using a simple Bayesian model, but it's not clear to me how to put that two standard deviations information into the model or into the standard error in such a way that I arrive at the right parameters in the end.

data {
  int<lower=0> J;
  real y[J];
  real<lower=0> sigma[J];

parameters {
  real theta[J];

model {
  y ~ normal(theta, sigma);

if you feel that the mean of your observations should be around 1100 and that 2500 is about two standard deviations away from the mean, then I think this quite accurately reflects your prior knowledge (i.e. mu should be centered around 1100 and sigma should have expected value close to (2500-1100)/2 = 700).

Looking at your Stan code, it seems like you consider sigma a part of your data, not of your parameters. That is appropriate if you know exactly what your standard deviation is going to be, but doesn't sound like your case. It seems more appropriate to put sigma in the parameters block. Second, it seems like both sigma and theta are now real arrays of length J. This means that every observation has its own mean and sigma, which is probably not what you want to accomplish.

The simplest choice of prior is completely uninformative, but this is rarely a practical choice of prior. A step up from there would be to choose a conjugate prior (see: https://en.wikipedia.org/wiki/Conjugate_prior - this page includes quite a few examples of priors for the Normal distribution). In this case, if you're using Stan, then you can pick almost any prior that reasonably reflects your assumptions. For instance, you could use a Normal prior for the mean and a Exponential prior for the standard deviation, each chosen such that they reflect your prior knowledge. For example, you could choose the following model to start, and then adjust as you see fit after checking its inferences.

  int<lower=0> J;
  real y[J];

  real mu;
  real<lower=0> sigma;

  y ~ normal(mu, sigma);
  mu ~ normal(1100, 100);
  sigma ~ exponential(1/700);

As you see, the prior for mu was chosen quite vaguely, with a wide range of plausible values (high standard deviation). One possible way to improve on this is if you have more specific knowledge that you can use. Or, you can set this as another parameter to be estimated from the data (essentially creating a hierarchical model). Another way to potentially improve on this model is to take a half-Cauchy prior for sigma (see Stan manual 2.15.0 page 124-128 for more on the choice of priors, or look here: https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations). These tend do do quite well in practice when you don't have much prior information. The problem in this case is that the Cauchy distribution doesn't have a defined expected value, so you can't use it to reflect your knowledge that sigma is 'about 700'.

I hope this helps. Just remember, regardless of the model you're trying to set up, if you have a dataset with only a single observation, it will always be difficult to learn a lot from the data.


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.