Setting response not as a distribution but rather as an equality

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\}:

blockNoise_i \sim \mathcal{N}(0,blockSD)

and

y_{l_j,i} = dose_{d_i}(1+blockNoise_{l_j})

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.

Just compute the slopes as transformed data and then fit a normal distribution to them.

data {
int<lower=0> N;
int<lower=0> nBlocks;
array[N] int blocks;
array[N] real y;
array[N] real doses;
}
transformed data {
vector[nBlocks] blockNoise;
for (i in 1:N) {
blockNoise[block[i]] = y[i]/doses[i] - 1;
}
}
parameters {
real<lower=0> blockSD;
}
model {
blockSD ~ normal(0,1);
blockNoise ~ normal(0,blockSD);
}