Skip to content

Commit 7e00829

Browse files
Implemented randomised PIT for discrete data, and included some more precomputed values for faster confidence bands.
1 parent c9ce1d3 commit 7e00829

File tree

2 files changed

+24
-18
lines changed

2 files changed

+24
-18
lines changed

R/ppc-distributions.R

Lines changed: 24 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -606,13 +606,13 @@ ppc_pit_ecdf <- function(y,
606606
if (is.null(pit)) {
607607
pit <- ppc_data(y, yrep) %>%
608608
group_by(.data$y_id) %>%
609-
dplyr::group_map(~ mean(.x$value[.x$is_y] >= .x$value[!.x$is_y])) %>%
609+
dplyr::group_map(
610+
~ mean(.x$value[.x$is_y] > .x$value[!.x$is_y]) +
611+
runif(1, max = mean(.x$value[.x$is_y] == .x$value[!.x$is_y]))
612+
) %>%
610613
unlist()
611614
if (is.null(K)) {
612-
K <- min(
613-
length(unique(ppc_data(y, yrep)$rep_id)) + 1,
614-
length(pit)
615-
)
615+
K <- nrow(yrep) + 1
616616
}
617617
} else {
618618
inform("'pit' specified so ignoring 'y', and 'yrep' if specified.")
@@ -632,20 +632,22 @@ ppc_pit_ecdf <- function(y,
632632
ggplot() +
633633
aes(
634634
x = 1:K / K,
635-
y = ecdf(pit)(seq(0, 1, length.out = K)) - (plot_diff == TRUE) * 1:K / K,
635+
y = ecdf(pit)(seq(0, 1, length.out = K)) -
636+
(plot_diff == TRUE) * seq(0, 1, length.out = K),
636637
color = "y"
637638
) +
638639
geom_step(show.legend = FALSE) +
639640
geom_step(aes(
640-
y = lims$upper[-1] / N - (plot_diff == TRUE) * 1:K / K,
641+
y = lims$upper[-1] / N - (plot_diff == TRUE) * seq(0, 1, length.out = K),
641642
color = "yrep"
642-
), show.legend = FALSE) +
643+
),
644+
linetype = 2, show.legend = FALSE) +
643645
geom_step(aes(
644-
y = lims$lower[-1] / N - (plot_diff == TRUE) * 1:K / K,
646+
y = lims$lower[-1] / N - (plot_diff == TRUE) * seq(0, 1, length.out = K),
645647
color = "yrep"
646-
), show.legend = FALSE) +
647-
yaxis_title(FALSE) +
648-
xaxis_title(FALSE) +
648+
),
649+
linetype = 2, show.legend = FALSE) +
650+
labs(y = "ECDF", x = "PIT") +
649651
yaxis_ticks(FALSE) +
650652
scale_color_ppc() +
651653
bayesplot_theme_get()
@@ -671,10 +673,13 @@ ppc_pit_ecdf_grouped <-
671673
if (is.null(pit)) {
672674
pit <- ppc_data(y, yrep, group) %>%
673675
group_by(.data$y_id) %>%
674-
dplyr::group_map(~ mean(.x$value[.x$is_y] >= .x$value[!.x$is_y])) %>%
676+
dplyr::group_map(
677+
~ mean(.x$value[.x$is_y] > .x$value[!.x$is_y]) +
678+
runif(1, max = mean(.x$value[.x$is_y] == .x$value[!.x$is_y]))
679+
) %>%
675680
unlist()
676681
if (is.null(K)) {
677-
K <- length(unique(ppc_data(y, yrep)$rep_id)) + 1
682+
K <- nrow(yrep) + 1
678683
}
679684
} else {
680685
inform("'pit' specified so ignoring 'y' and 'yrep' if specified.")
@@ -723,13 +728,14 @@ ppc_pit_ecdf_grouped <-
723728
geom_step(aes(
724729
y = .data$lims_upper - (plot_diff == TRUE) * .data$x,
725730
color = "yrep"
726-
), show.legend = FALSE) +
731+
),
732+
linetype = 2, show.legend = FALSE) +
727733
geom_step(aes(
728734
y = .data$lims_lower - (plot_diff == TRUE) * .data$x,
729735
color = "yrep"
730-
), show.legend = FALSE) +
731-
xaxis_title(FALSE) +
732-
yaxis_title(FALSE) +
736+
),
737+
linetype = 2, show.legend = FALSE) +
738+
labs(y = "ECDF", x = "PIT") +
733739
yaxis_ticks(FALSE) +
734740
bayesplot_theme_get() +
735741
facet_wrap("group") +

R/sysdata.rda

991 Bytes
Binary file not shown.

0 commit comments

Comments
 (0)