I'm taking a class and all of the code/modeling examples are done in WinBugs. However, the professor has said that we can use whichever library/language we want when we turn in the homework.
I'm trying to convert the WinBugs examples from class into Stan, but I'm having problems getting my Stan model file correct.
model{
for(i in 1:N){
y[i] ~ dnorm(mu[i], tau)
mu[i] <- beta0 + beta1*x[i]
}
}
#priors
beta0 ~ dnorm(0,0.001)
beta1 ~ dnorm(0,0.001)
tau ~ dgamma(0.01, 0.01)
#prediction
new.x <- 90
ympred <- beta0 + beta1 * new.x
ypred ~ dnorm(ympred, tau)
I've come up with the following code for my stan file.
data {
int<lower=0> N;
vector[N] dependent_variable_y;
vector[N] independent_variable_x;
}
// these are
parameters {
real beta0;
real beta1;
real tau;
real mu;
}
// The model to be estimated. We model the output
// 'y' to be normally distributed with mean 'mu'
// and standard deviation 'sigma'.
model {
//priors
beta0 ~ normal(0, 0.001);
beta1 ~ normal(0, 0.001);
tau ~ gamma(0.01, 0.01);
for (i in 1:N){
mu = beta0 + beta1 * independent_variable_x[i]
}
}
The error is currently telling me that I'm missing a } somewhere, but I can't figure out where.
You're missing a semicolon in the last relevant line
mu = beta0 + beta1 * independent_variable_x[i];
That's probably where the missing } is coming from. But more importantly, there's a few mistakes in the model. First, mu is determined completely by the beta
s and x
, it's not an independent parameter. Either you put it in the transformed parameters
block, or just omit saving it altogether. Second, there's no error model. You're not defining how y
depends on mu
and sigma
(assuming tau
is supposed to be sigma?). Third, (not a mistake, just efficiency) you don't have to loop over N
, you can just vectorize. Fourth, you're using EXTREMELY narrow priors on the betas, those priors pretty much say that you are sure the betas are 0. You should have a good reason to do so. Assuming you do have that, here's a version that should work:
data {
int<lower=0> N;
vector[N] y;
vector[N] x;
}
parameters {
real beta0;
real beta1;
real tau;
}
model {
//priors
beta0 ~ normal(0, 0.001);
beta1 ~ normal(0, 0.001);
tau ~ gamma(0.01, 0.01);
y ~ normal(beta0 + beta1 * x, tau);
}