I am trying to fit a probit model where each observation happened on a specific *dose* \{x_{1} \ldots x_m\} and for a specific block \{b_1 \ldots b_K\}. We also assume that each block has block-specific noise, \{blockNoise_1 \ldots blockNoise_K\} that affects the slope, where:

The probit model would be simulated as follows:

Where \Phi is the normal cumulative distribution function.

Given the \{y_{i,j}\}, I am trying to infer \text{block_SD}.

However, I have been really struggling to get this to work, things work fine computationally but I usually only get the mean of the prior **0.1** as an answer:

```
data {
int<lower=0> N;
int<lower=0> nBlocks;
array[N] int blocks;
array[N] int<lower=0, upper=1> y;
array[N] real doses;
}
parameters {
real<lower=0.01,upper=0.21> blockSD;
vector[nBlocks] blockSDs;
}
transformed parameters {
array[N] real probits;
vector[N] p;
for (i in 1:N) {
probits[i] = doses[i]*(blockSDs[blocks[i]] + 1);
p[i] = Phi(probits[i]);
}
}
model {
blockSD ~ normal(0.1,0.05);
blockSDs ~ normal(0,blockSD);
y ~ bernoulli(p);
}
```

(for the model to work reliably the parameters need to be non-centered but I kept it like this for simplicity)

I have wondered if itâ€™s the data that I am generating - in that there is just not enough information to infer anything. Below is an example (although I use more replicates and blocks and the result is similar). In the figure below, on top is the blockSD=0.001 (i.e the slope in the model does not change much between blocks), but it is higher on the bottom (0.2):

and I also plot the probits:

(obviously there is difference in the variation in the slopes).

What am I missing? I think the problem is in the way that I am specifying the model but I am not sure what the issue is, it seems like itâ€™s not sampling things correctly.