# Linear model with generated quantities - ypred issues?

**URL:** <https://discourse.mc-stan.org/t/linear-model-with-generated-quantities-ypred-issues/2239>\
**Category:** Modeling\
**Tags:** fitting-issues\
**Created:** [October 18, 2017, 4:56pm UTC](https://discourse.mc-stan.org/t/linear-model-with-generated-quantities-ypred-issues/2239 "2017-10-18T16:56:28Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ara\_Winter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ara_winter/32/4284_2.png) [@Ara\_Winter](https://discourse.mc-stan.org/u/Ara_Winter)\
**Post date:** [October 18, 2017, 4:56pm UTC](https://discourse.mc-stan.org/t/linear-model-with-generated-quantities-ypred-issues/2239/1 "2017-10-18T16:56:28Z")

</div>

Hello all,

I am still slowly moving along learning PyStan and shifting our models over. I am trying to understand how the generated quantities work. I have a simple linear model of ocean fish catch plotted against ocean net primary production (NPP). Previous work shows that this can be modeled roughly with a simple linear model. I also want to see how well the model does. I added a generated quantities section to the model.

I log transformed the data. Which thinking about this might be my issue with the model.

When I run fit.plot([‘y\_pred’]) that values look weird to me. They range 4 to 6 while the original data goes from -8 to 2.

The juypter notebook and raw data are [here](https://github.com/bioinfonm/pystan_musings/tree/master/beginners_excercise/fish_catch_productivity)

The model I used was this:

> regress\_code = “”"
> 
> data {  
> int\<lower=0\> N;  
> int N2;  
> vector[N] y;  
> vector[N] x;  
> }
> 
> parameters {  
> real alpha;  
> real beta;  
> real\<lower=0\> sigma;  
> }
> 
> model {  
> alpha ~ normal(0,10);  
> beta ~ cauchy(0,5);  
> sigma ~ cauchy(0,5);  
> for (n in 1:N)  
> y[n] ~ normal(alpha + beta \* x[n], sigma);  
> }
> 
> generated quantities {  
> real y\_pred;
> 
> y\_pred = normal\_rng(alpha + beta \* x[N2], sigma);
> 
> }
> 
> “”"

I picked cauchy for beta since in the worse case the slope could be zero but likely it’s some positive number.

Did I set my priors up correctly and ypred? Or should y\_pred not be normal\_rng? Should I log transform before the model or specify that in the model?

Thank you!  
ara

---

<div class="post-metadata">

**Author:** ![Ara\_Winter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ara_winter/32/4284_2.png) [@Ara\_Winter](https://discourse.mc-stan.org/u/Ara_Winter)\
**Post date:** [October 18, 2017, 9:08pm UTC](https://discourse.mc-stan.org/t/linear-model-with-generated-quantities-ypred-issues/2239/2 "2017-10-18T21:08:10Z")

</div>

I may have made some headway after diving back into the Stan manual and some old Stan forum posts.

---

<div class="post-metadata">

**Author:** ![HeltonMaciel](https://avatars.discourse-cdn.com/v4/letter/h/5f9b8f/32.png) [@HeltonMaciel](https://discourse.mc-stan.org/u/HeltonMaciel)\
**Post date:** [October 19, 2017, 9:47am UTC](https://discourse.mc-stan.org/t/linear-model-with-generated-quantities-ypred-issues/2239/3 "2017-10-19T09:47:46Z")

</div>

> [@Ara\_Winter](#):
>
> > generated quantities {  
> > real y\_pred;
> > 
> > y\_pred = normal\_rng(alpha + beta \* x[N2], sigma);
> > 
> > }

looks like you’re using only one value of the covariate `x` by setting `x[N2]` outside a loop.

i’d do something like this:

> [@Ara\_Winter](#):
>
> > regress\_code = “”"
> > 
> > data {  
> > int\<lower=0\> N;  
> > int N2;  
> > vector[N] y;  
> > vector[N] x;  
> > }
> > 
> > parameters {  
> > real alpha;  
> > real beta;  
> > real\<lower=0\> sigma;  
> > }
> > 
> > transformed parameters{  
> > vector[N] mu
> > 
> > for (n in 1:N){  
> > mu[n] = alpha + beta\*x[n]  
> > }  
> > }  
> > model {  
> > alpha ~ normal(0,10);  
> > beta ~ cauchy(0,5);  
> > sigma ~ cauchy(0,5);
> > 
> > ```
> > y ~ normal(mu, sigma);
> > 
> > ```
> > 
> > }
> > 
> > generated quantities {  
> > vector[N] y\_pred;
> > 
> > for (n in 1:N){  
> > y\_pred[n] = normal\_rng(mu[n], sigma);  
> > }
> > 
> > }
> > 
> > “”"

---

<div class="post-metadata">

**Author:** ![Ara\_Winter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ara_winter/32/4284_2.png) [@Ara\_Winter](https://discourse.mc-stan.org/u/Ara_Winter)\
**Post date:** [October 19, 2017, 12:47pm UTC](https://discourse.mc-stan.org/t/linear-model-with-generated-quantities-ypred-issues/2239/4 "2017-10-19T12:47:51Z")

</div>

Ok, let me try this and see if I can figure out how all this is working.

Thanks for taking the time!

Ara

---

<div class="post-metadata">

**Author:** ![Ara\_Winter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ara_winter/32/4284_2.png) [@Ara\_Winter](https://discourse.mc-stan.org/u/Ara_Winter)\
**Post date:** [October 27, 2017, 9:19pm UTC](https://discourse.mc-stan.org/t/linear-model-with-generated-quantities-ypred-issues/2239/5 "2017-10-27T21:19:25Z")

</div>

@HeltonMaciel Thanks! I think I have it working. If I understand correctly I was only generating one y or a given x. With the for loop I get all the y’s for a given x.

I also went back to just a simple linear model and then built up the model to a generated quantity model.

Ara

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [November 2, 2017, 10:36pm UTC](https://discourse.mc-stan.org/t/linear-model-with-generated-quantities-ypred-issues/2239/6 "2017-11-02T22:36:03Z")

</div>

Do the `y_pred` not look like the original `y` even if you simulate the original `y` according to the model?

---

<div class="post-metadata">

**Author:** ![Ara\_Winter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ara_winter/32/4284_2.png) [@Ara\_Winter](https://discourse.mc-stan.org/u/Ara_Winter)\
**Post date:** [November 2, 2017, 10:51pm UTC](https://discourse.mc-stan.org/t/linear-model-with-generated-quantities-ypred-issues/2239/7 "2017-11-02T22:51:56Z")

</div>

In the raw data y goes from 0.000305 to 9.3 (mean ten year fish catch, tons per km2).

y\_pred (mean) only covers part of this range 0.39 to 2.63. The max y\_pred covers the upper range but the min y\_pred goes negative. I found that weird since you can’t have a negative tonnage.

This leaves me wondering if I didn’t set up the model or priors correctly. Or the original data is just not that great.

My notebook for this is on github [here](https://github.com/bioinfonm/pystan_musings/blob/master/beginners_excercise/fish_catch_productivity/pystan_linear_model_fish_production_from_npp.ipynb).

Raw data is [here](https://github.com/bioinfonm/pystan_musings/blob/master/beginners_excercise/fish_catch_productivity/fish_catch_LME_npp_cb_per_km2.txt).

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [November 2, 2017, 10:55pm UTC](https://discourse.mc-stan.org/t/linear-model-with-generated-quantities-ypred-issues/2239/8 "2017-11-02T22:55:17Z")

</div>

That’s why I was suggesting simulating data from the model, then fitting it. Then the `y_pred` should look right. It’s the only way to separate model misspecification from implementation issues.

---

<div class="post-metadata">

**Author:** ![Ara\_Winter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ara_winter/32/4284_2.png) [@Ara\_Winter](https://discourse.mc-stan.org/u/Ara_Winter)\
**Post date:** [November 2, 2017, 11:00pm UTC](https://discourse.mc-stan.org/t/linear-model-with-generated-quantities-ypred-issues/2239/9 "2017-11-02T23:00:44Z")

</div>

Oh, sorry. I misunderstood. Let me try that.

Thanks!  
ara

---

<div class="post-metadata">

**Author:** ![Ara\_Winter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ara_winter/32/4284_2.png) [@Ara\_Winter](https://discourse.mc-stan.org/u/Ara_Winter)\
**Post date:** [November 3, 2017, 1:21pm UTC](https://discourse.mc-stan.org/t/linear-model-with-generated-quantities-ypred-issues/2239/10 "2017-11-03T13:21:28Z")

</div>

Just to make sure I am doing this right. I hit up the Stan manual 2.16 and on page 160 under Posterior Predictive Checks.

I set up in my model:

generated quantities {  
vector[N] y\_tilde;  
for (n in 1:N)  
y\_tilde[n] = normal\_rng(alpha + beta \* x[n], sigma);  
}

EDIT: Answered my own question about the tilde.

---

<div class="post-metadata">

**Author:** ![Ara\_Winter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ara_winter/32/4284_2.png) [@Ara\_Winter](https://discourse.mc-stan.org/u/Ara_Winter)\
**Post date:** [November 3, 2017, 1:40pm UTC](https://discourse.mc-stan.org/t/linear-model-with-generated-quantities-ypred-issues/2239/11 "2017-11-03T13:40:06Z")

</div>

So I ended up with:

Random data

```
x_tilde = random.sample(range(50, 350), 100)
N_tilde = len(x_tilde)

np.mean(x_tilde, axis=0)
192.27000000000001

```

Dictionary for Stan

```
Make a dictionary containing all data to be passed to STAN
regress_sim_dat = {'x': x,
 'y': y,
 'N': N,
 'x_tilde':x_tilde,
 'N_tilde':N_tilde}

```

The model, priors, and generated data.

```
 regress_sim = """
 data {
   int<lower=0> N;
   vector[N] y;
   vector[N] x;
   int<lower=0> N_tilde;
   vector[N_tilde] x_tilde;
 }
 parameters {
   real alpha;
   real beta;
   real<lower=0> sigma;
 }
 model {
   alpha ~ normal(0,10);    
   beta ~ cauchy(0,5);
   sigma ~ cauchy(0,5);
 
   for (n in 1:N)
     y[n] ~ normal(alpha + beta * x[n], sigma);
 }
 generated quantities {
    vector[N_tilde] y_tilde;
       
    for (n in 1:N_tilde)
      y_tilde[n] = normal_rng(alpha + beta * x_tilde[n], sigma);
 }
 """

```

EDIT: Cleaned up the code blocks with ```

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [November 3, 2017, 1:55pm UTC](https://discourse.mc-stan.org/t/linear-model-with-generated-quantities-ypred-issues/2239/12 "2017-11-03T13:55:55Z")

</div>

Yes, that looks right. To set off code, you want triple back ticks on their own line before and after:

```
```
your code
```
```

---

<div class="post-metadata">

**Author:** ![Ara\_Winter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ara_winter/32/4284_2.png) [@Ara\_Winter](https://discourse.mc-stan.org/u/Ara_Winter)\
**Post date:** [November 3, 2017, 3:33pm UTC](https://discourse.mc-stan.org/t/linear-model-with-generated-quantities-ypred-issues/2239/13 "2017-11-03T15:33:20Z")

</div>

Ok, cool. I think I understand all this. The y\_tilde range with the model above gives a y\_tilde mean between ~0.4 to ~2.8 . There are still negative numbers in the y\_tilde traceplot.

 ![fish_npp_y_tilde_sim_data](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/0/04f121047a339e56eb505065d1b393410711d7aa.png)