Home Ask Login Register

Developers Planet

Your answer is one click away!

spore234 February 2016

simple Gamma GLM in STAN

I'm trying a simple Gamma GLM in STAN and R, but it crashes immediately

generate data:

N<-500 #sample size
#the model
K<-dim(X)[2] #number of regression params
#the regression slopes
shape <- 10
#simulate gamma data

this is my STAN model:

model_string <- "
data {
  int<lower=0> N; //the number of observations
  int<lower=0> K; //the number of columns in the model matrix
  matrix[N,K] X; //the model matrix
  vector[N] y; //the response
parameters {
  vector[K] betas; //the regression parameters
  real<lower=0, upper=1000> shape; //the shape parameter
model {  
  y ~ gamma(shape, (shape/exp(X * betas)));

when I run this model, R immediately crashes:

m <- stan(model_code = model_string, data = list(X=X, K=3, N=500, y=y), chains = 1, cores=1)

update : I think the problem is somewhere in the vectorization as I can get a running model where I pass every column of X as a vector.

update2: this also works

  for(i in 1:N)
    y[i] ~ gamma(shape, (shape / exp(X[i,] * betas)));


Ben Goodrich February 2016

The problem with the original code is that there is no operator currently defined in Stan for a scalar divided by a vector. In this case, shape / exp(X * betas) You might be able to do shape[1:N] ./ exp(X * betas) or failing that, (shape * ones_vector) ./ exp(X * betas)

Post Status

Asked in February 2016
Viewed 2,446 times
Voted 4
Answered 1 times


Leave an answer

Quote of the day: live life