# Help implementing a Gamma model (rstan)

**URL:** https://discourse.mc-stan.org/t/help-implementing-a-gamma-model-rstan/37475
**Category:** Modeling
**Tags:** specification, fitting-issues, rstan
**Created:** [December 2, 2024, 3:53pm UTC](https://discourse.mc-stan.org/t/help-implementing-a-gamma-model-rstan/37475 "2024-12-02T15:53:28Z")
**Posts on this page:** 3
**Page:** 1

<div class="post-metadata">

### Author: ![Ruben222](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ruben222/32/5514_2.png) [@Ruben222](https://discourse.mc-stan.org/u/Ruben222)
#### Post date: [December 2, 2024, 3:53pm UTC](https://discourse.mc-stan.org/t/help-implementing-a-gamma-model-rstan/37475/1 "2024-12-02T15:53:29Z")

</div>

Hello. I am having a lot of trouble trying to get a Gamma model working with some data.

A while ago, I wrote myself a simple tutorial model to learn how gamma models worked (shown below).

```stan
data {
  int<lower=1> n; // Number of observations
  vector[n] y; // Response variable
  vector[n] x; // Predictor variable
}

parameters {
  real a;         
  real b;                
  real<lower=0> shape;    
}

transformed parameters {
  vector[n] mu; 
  vector[n] beta; /
  for(i in 1:n){
    mu[i] = exp(a + b * x[i]); 
    beta[i] = shape / mu[i]; . 
  }
}

model {
  // Likelihood
  y ~ gamma(shape, beta); 
  // Priors
  a ~ normal(0, 2);
  b ~ normal(0, 2);
  shape ~ gamma(0.01, 0.01);
}

```

This model correctly recovers values from fake data, so I think it is working as intended.

```no-highlight
N <- 100
x <- runif(N, -1, 1)
a <- 0.5
b <- 1.1
y_true <- exp(a + b * x)
shape <- 10
y <- rgamma(N, rate = shape / y_true, shape = shape)

fake_data_a <- data.frame(x=x,y=y)

```

So, I’ve been trying to run a model specified in the same way with some real data. I.e.

```stan
data {
  int<lower=1> n;       
  vector[n] G_value;            
}

parameters {
  real a;
  real<lower=0> shape;  
}

transformed parameters {
  vector[n] mu; 
  vector[n] beta; 
  for(i in 1:n){
    mu[i] = exp(a); 
    beta[i] = shape / mu[i]; 
  }
}

model {
  G_value ~ gamma(shape, beta); 
  a ~ normal(0,5);
  shape ~ gamma(0.01, 0.01);
}

```

However, no matter what I try the model never samples. Irrespective of the predictors that are included. Above I simplified it as much as possible. It always produces errors like this:

```no-highlight
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can't start sampling from this initial value

```

I’m sure it will be something quite obvious, but I can’t figure out what is going wrong. Does anybody have any suggestions?

For example, is there anything wrong with how I have setup either of my models?

Thanks in advance.

---

<div class="post-metadata">

### Author: ![Garren\_Hermanus](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/garren_hermanus/32/17356_2.png) [@Garren\_Hermanus](https://discourse.mc-stan.org/u/Garren_Hermanus)
#### Post date: [December 2, 2024, 4:21pm UTC](https://discourse.mc-stan.org/t/help-implementing-a-gamma-model-rstan/37475/2 "2024-12-02T16:21:46Z")

</div>

Hi. This is a simple model, hence stan shouldn’t have any issues with initialization. I would check that the data you are feeding is actually positive.

> [@Ruben222](#):
>
> ```no-highlight
> data {
> int<lower=1> n;       
> vector[n] G_value;            
> }
> 
> parameters {
> real a;
> real<lower=0> shape;  
> }
> 
> transformed parameters {
> vector[n] mu; 
> vector[n] beta; 
> for(i in 1:n){
> mu[i] = exp(a); 
> beta[i] = shape / mu[i]; 
> }
> }
> 
> model {
> G_value ~ gamma(shape, beta); 
> a ~ normal(0,5);
> shape ~ gamma(0.01, 0.01);
> }
> 
> ```

Also, for performance, most of the stan functions are optimized for vectorization. If data follows the same distribution, you can define the data as arrays and feed scalar values for the parameters. The optimized code, with the addition of checking for positive data would be

```no-highlight
data {
  int<lower=1> n;       
  array[n] real<lower=0> G_value;            
}

parameters {
  real a;
  real<lower=0> shape;  
}

transformed parameters {
  real mu = exp(a); 
  real beta = shape / mu; 
}

model {
  G_value ~ gamma(shape, beta); 
  a ~ normal(0,5);
  shape ~ gamma(0.01, 0.01);
}

```

---

<div class="post-metadata">

### Author: ![Ruben222](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ruben222/32/5514_2.png) [@Ruben222](https://discourse.mc-stan.org/u/Ruben222)
#### Post date: [December 4, 2024, 8:39am UTC](https://discourse.mc-stan.org/t/help-implementing-a-gamma-model-rstan/37475/3 "2024-12-04T08:39:00Z")

</div>

Thank you very much for your reply. You were right, there were a few 0s in the data which was stopping the model running. Removing them let everything work as intended.

Thanks also for posting a vectorised version of the model. I’ve been able to expand it to include predictors without needing the for loop, which has sped things up a lot.