This is probably very simple but I have been thinking about it too long and I don’t know where to look. Have some experience with Stan but have never come across this before

I am playing with trying to fit a simple model where I generate a set of lines \{l_1 \ldots l_K\}, where the slope of each line differs by a corresponding block-specific noise, \{blockNoise_1 \ldots blockNoise_K\}:

and

where *doses* are dose_{1} \ldots dose_m

I can easily generate data as:

```
data {
int<lower=0> N;
int<lower=0> nBlocks;
array[N] real doses;
array[N] int blocks;
}
generated quantities {
real<lower=0> blockSD;
vector[nBlocks] blockNoise;
vector[N] y;
blockSD = abs(normal_rng(0,0.1));
for (i in 1:nBlocks) {
blockNoise[i] = normal_rng(0,blockSD);
}
for (i in 1:N) {
y[i] = doses[i]*(1+blockNoise[blocks[i]]);
}
}
```

The question is, how do I fit this? My goal is to estimate blockSD. I am not sure how to use “y”, or maybe more generally, if the response does not exactly follow a distribution but is rather a function of elements that follow a distribution how does one use it in the program?

```
data {
int<lower=0> N;
int<lower=0> nBlocks;
array[N] int blocks;
array[N] real y;
array[N] real doses;
}
parameters {
real<lower=0> blockSD;
vector[nBlocks] blockNoise;
}
model {
blockSD ~ normal(0,1);
blockNoise ~ normal(0,blockSD_b);
for (i in 1:N) {
y[i] = doses[i]*(1+blockNoise[blocks[i]]); // THIS DOES NOT WORK!
}
}
```

Any help appreciated.