LinearModel.stan (817 Bytes) LinearModel_stan.R (1.3 KB)
Hi,
I’m trying to implement a linear model with Zellner prior:
Y | X,\beta, \sigma^2 \sim N(X \beta, \sigma^2 I_n)
\beta | \sigma^2, X \sim N_p(b_0, \sigma^2 B_0)
\sigma^2 \sim Inv-gamma(\frac{\nu0}{2}, \frac{\nu_0\sigma_0^2}{2})
with \nu_0=0 and B_0=c (X^T X)^{-1} (c constant)
To run my model I used a simulated dataset where
\beta=[1,2,3,4]
and
\sigma^2=5
The \beta estimated by the model seems to converge to the true value (only the intercept shows a slightly different value), but the variance doesn’t converge to the exact value.
This is my model
data {
int<lower=0> N; // number of data
vector[N] y; // response
int<lower=0> p; //number of regressors
matrix[N, p] X; //matrix of covariates
// Zellner prior parameters
real c;
vector[p] b0;
// prior parameters of gamma
real sigma0;
real nu0;
}
parameters {
vector[p] beta; //regressors paramters (column vector)
real<lower=0> sigma2; // standard deviation
}
transformed parameters
{
vector[N] mu; // mean
matrix[p, p] B0;
for(i in 1:N){
mu[i] = X[i,] * beta;
}
B0= c * inverse(X' * X);
}
model {
//likelihood:
y ~ normal(mu, pow(sigma2, 0.5));
//Prior:
beta ~ multi_normal(b0, B0);
sigma2 ~ inv_gamma(nu0/2, nu0*sigma0/2); // paramters of zellner prior
}
}
This is what I get
Inference for Stan model: LinearModel.
4 chains, each with iter=1e+05; warmup=50000; thin=10;
post-warmup draws per chain=5000, total post-warmup draws=20000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
beta[1] 1.42 0.02 3.34 -5.07 -0.83 1.43 3.69 7.97 19742 1
beta[2] 2.02 0.00 0.28 1.48 1.83 2.02 2.20 2.56 19823 1
beta[3] 2.92 0.00 0.17 2.59 2.81 2.92 3.03 3.26 19299 1
beta[4] 3.95 0.00 0.15 3.66 3.85 3.95 4.05 4.24 19493 1
sigma2 **112.68** 0.12 16.3 85.28 101.11 111.35 122.68 148.58 19549 1
Samples were drawn using NUTS(diag_e) at Thu Dec 19 18:01:11 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
I also tried to use reference priors and i get an estimated value of \sigma^2 equal to 6.34 which is more closed to what I expect.
Can someone help me to understand why using the zellner prior I get that result? Is there something wrong with my model?
Thank you in advance!