Search code examples
machine-learningstatisticsbayesianstanwinbugs

Am I defining my model parameters correctly when converting from WinBugs to Stan?


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.


Solution

  • 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 betas 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);
            
    }