Stan model for noisy observations of a sigmoidal function, where the noise is truncated to a finite interval

TruncSigmoid.stan

functions {
  vector ft(vector t, real r, real y0, real tC, real L) {
    vector[num_elements(t)] mu;
    vector[num_elements(t)] exprt;
    exprt = exp(r*(t-tC));
    mu = y0*L*exprt ./ (L + y0*(exprt - 1));
    return mu;
  }
  
  vector dfdt(vector t, real r, real tC, real L) {
    vector[num_elements(t)] dmu;
    vector[num_elements(t)] sddenom;
    for (i in 1:num_elements(t)) {
      sddenom[i] = ((exp(r*tC) + exp(r*t[i]))^(-2));
    }
    dmu = r*L*exp(r*(t+tC)) .* sddenom;
    return dmu;
  }
}
data {
  int<lower = 1> M;
  int<lower = 1> N;
  real<lower = 1> maxY;
  real tcrit;
  matrix<lower=0, upper=maxY>[M,N] y;
  vector[M] t;
  real eps;
  int<lower = 1> P;
  vector[P] tpred;
}
parameters {
  real<lower = 0> r;
  real<lower = 0, upper=maxY> mu0;
}
transformed parameters {
  vector[M] curr_mu;
  vector[M] curr_sd;
  curr_mu = ft(t, r, mu0, tcrit, maxY);
  curr_sd = dfdt(t, r, tcrit, maxY) + eps;
}
model {
  r ~ gamma(1, 10);
  mu0 ~ uniform(0, maxY);
  for (i in 1:M) {
    y[i,] ~ normal(curr_mu[i], curr_sd[i]);
  }
}
generated quantities{
  vector[P] pred_mu;
  vector[P] pred_sd;
  pred_mu = ft(tpred, r, mu0, tcrit, maxY);
  pred_sd = dfdt(tpred, r, tcrit, maxY) + eps;
}