I have a 2-step generative model.

- vector i ~ normal(mu, std)
- I algebraically transform each of the entries of vector i onto a scale of 0-1 using a sigmoid function.
- Then I draw a new variable vector y ~ Bernoulli( sigmoid(i) )
- The observed data, vector s, is equal to: { i if y=1 } and { 0 if y=0 }

I hope to draw samples for the parameters mu and std, given I observe the vector of data s. I have coded the model as follows, but it is giving me a lot of divergences. (Even though I am following non-centered (centered?) implementation.) I am also not sure if I have coded my generative process correctly.

```
"""
functions {
real sigm(real x){
return (1 / (1 + exp(-0.92*(x-8.24)));
}
real signal_lpdf(real s, real i, real mu, real std){
if (s==0) {
return bernoulli_lpmf(0 | sigm(i)) + normal_lpdf((i-mu)/std | 0, 1);
}
else {
return bernoulli_lpmf(1 | sigm(s))+ normal_lpdf((s-mu)/std | 0, 1);
}
}
}
data {
int<lower=0> nPeps; // number of theoretical peptides in protein
real<lower=0> s[nPeps]; // array holding signal strength for that protein (0 is missing value)
}
parameters {
real mu;
real i;
real<lower=0> std;
}
transformed parameters {
real mu2;
real<lower=0> std2;
mu2 = mu*10+7.15;
std2 = std*10+3;
}
model {
mu ~ std_normal();
std ~ std_normal();
for (j in 1:nPeps) {
s[j] ~ signal(i, mu2, std2);
}
}
"""
```