Specifying a conditional prior

Modeling
, ,
1

I haven’t been able to find any help files/manuals on this topic, but wanted to see if it is possible to specify a prior in STAN (or brms) such that part of the prior is a conditional distribution? For example, in the case of using a normal inverse gamma in regression, you may want to decompose the prior for beta and sigma^2 as

p(beta, sigma^2) = p(beta|sigma^2)p(sigma^2)

Is there a way to specify the prior to capture that conditional distribution?

To make the example somewhat more concrete, suppose I have data as follows:

n <- 10
beta0 <- 5
beta1 <- 3.7
x <- 1:n
y <- beta0 + beta1 * x + rnorm(n)
dat <- data.frame(y,x)

X <- cbind(1, x)
betahat <- solve(t(X) %*% X) %*% t(X) %*% y

So my statistical model is Y = betaX + epsilon where beta = (beta0, beta1) and epsilon ~ N(0, sigma^2) (so the classical regression setup). I need to place priors on beta and sigma^2, and want my prior to be of the form

p(beta, sigma^2) = p(beta|sigma^2) p(sigma) where

beta | sigma^2 ~ N( betahat, sigma^2 * (X’X)^{-1})
sigma^2 ~ Inverse Gamma(5, 10)

2

Sure! For example:

parameters{
 real x;
 real y;
}
model{
 x ~ std_normal();
 y ~ normal(x, 1);
}

The above example makes the conditional structure explicit, but it would be just as valid to figure out the joint prior and to write it down in one statement (in this case it’d be multivariate normal).

In your example, beware the need for a Jacobian adjustment if you declare parameter sigma and put a prior on sigma^2.

3

@jsocolar Thanks for the prompt reply. I updated my question with some actual data to be more concrete. I wonder if you could put your solution in the context of my multivariate normal and inverse gamma priors?

4

I’m having trouble understanding what parts of the data simulation are supposed to also be encoded in the model, and what parts are just a way of simulating some data. In the model, is betahat a transformed parameter or is it passed as data?

If betahat is a transformed parameter, and the transformation mirrors the way that it’s computed in the simulate data, then you have that betahat depends on y, and y depends on beta, but that the prior for beta depends on betahat. Seems like some heavy math lifting would be required to turn this into a generative prior that’s easy to reason about.

If betahat is data in your model, then:

model{
 sigma2 ~ inverse_gamma(5, 10);
 beta ~ normal(betahat, sigma2*inverse(X' * X));
}

Again, remember that a Jacobian adjustment will be required if you choose to declare parameter sigma and derive the variance as a transformation of the parameter, rather than just directly parameterizing in terms of the variance as above.

5

@jsocolar betahat is being passed as data.

Here is an update with R and Stan code:

library(rstan)

n <- 100
beta0 <- 5
beta1 <- 3.7
x <- 1:n
y <- beta0 + beta1 * x + rnorm(n, 0, 2.4)
dat <- data.frame(y,x)

X <- cbind(1, x)
k <- ncol(X)
betahat <- solve(t(X) %*% X) %*% t(X) %*% y
s2 <- c(t(y - X %*% betahat) %*% (y - X %*% betahat) / (n - k))

my_data <- list(
  n = n,
  y = y,
  X = X,
  k = k,
  mu0 = c(betahat),
  delta0 = solve(t(X) %*% X),
  a0 = (n - k) / 2,
  b0 = (n - k) * s2 / 2
) 

fit <- stan(
  file = "linear model.stan",
  data = my_data,
  chains = 1,
  warmup = 10000,
  iter = 100000
)

Stan Code:

data {
  // the input data
  int<lower = 1> n;               // total number of obs
  vector[n] y;                    // response vector
  int<lower = 1> k;               // number of coefficients
  matrix[n, k] X;                 // design matrix

  // parameters for the priors
  vector[k] mu0;
  matrix[k, k] delta0;
  real<lower = 0> a0;
  real<lower = 0> b0;
}

parameters {
  vector[k] beta;
  real<lower = 0> sigma2;
}

model {
 // prior
 sigma2 ~ inv_gamma(a0, b0);
 beta ~ multi_normal(mu0, sigma2 * delta0);

 // likelihood
 y ~ normal(X * beta, sigma2);
 
 }

I think this works now, the only thing I am unsure of is if I am passing the right terms to the likelihood. sigma2 is a variance with an inverse gamma prior, but do I have to pass a standard deviation to the y ~ normal() command?

6

Yes, pass a standard deviation to the y ~ normal! Sorry that wasn’t clear.

7

@jsocolar is it as simple as just updating the code to the following:

y ~ normal(X * beta, sqrt(sigma2));

Likewise, does multi_normal also expect a standard deviation for the covariance matrix? As in,

beta ~ multi_normal(mu0, sigma2 * delta0);

Is passing a variance multiplied by a matrix.