joint_reporting_delay_with_r_onsets <- joint_rt_fit |>
gather_draws(reporting_delay[day]) |>
ungroup() |>
filter(.draw %in% sample(.draw, 100)) |>
mutate(day = day - 1) # Stan's indexing starts from 1, not zero!
reporting_delay_pmf_df = cbind(day=0:15, pmf=reporting_delay_pmf)
ggplot() +
geom_col(
data=reporting_delay_pmf_df, mapping=aes(x = day, y = reporting_delay_pmf), alpha = 0.6
) +
geom_line(
data=joint_reporting_delay_with_r_onsets, mapping = aes(x = day, y = .value, group = .draw), alpha = 0.1, color="red") +
labs(title = "Reporting delay distribution, estimated (red) and observed (gray)", y="PMF")
