stancode=
"
data{
real onset_delay;
real ln_par1;
real ln_par2;
}
parameters{
real <lower=0> d_incubation;
real <lower=0> c_incubation;
}
transformed parameters{
real infection_delay = onset_delay + d_incubation - c_incubation;
}
model{
d_incubation ~ lognormal(ln_par1,ln_par2);
c_incubation ~ lognormal(ln_par1,ln_par2);
}

"
model=stan_model(model_code = stancode)
fn.stan <- function(delay){
dataList=list(onset_delay=delay,ln_par1=ln.par1,ln_par2=ln.par2)
fit=sampling(model,data=dataList)
ms=rstan::extract(fit)
mean(ms$infection_delay < 0)
}
fn.stan(2)
onset_delays=0:20
y=sapply(onset_delays,fn.stan)
plot(onset_delays,y, ylab='Pr[ Infected Later ])',axes=F) ; axis(1)