Skip to content

Commit

Permalink
Implemented randomised PIT for discrete data, and included some more …
Browse files Browse the repository at this point in the history
…precomputed values for faster confidence bands.
  • Loading branch information
TeemuSailynoja committed Jan 19, 2024
1 parent c9ce1d3 commit 7e00829
Show file tree
Hide file tree
Showing 2 changed files with 24 additions and 18 deletions.
42 changes: 24 additions & 18 deletions R/ppc-distributions.R
Original file line number Diff line number Diff line change
Expand Up @@ -606,13 +606,13 @@ ppc_pit_ecdf <- function(y,
if (is.null(pit)) {
pit <- ppc_data(y, yrep) %>%
group_by(.data$y_id) %>%
dplyr::group_map(~ mean(.x$value[.x$is_y] >= .x$value[!.x$is_y])) %>%
dplyr::group_map(
~ mean(.x$value[.x$is_y] > .x$value[!.x$is_y]) +
runif(1, max = mean(.x$value[.x$is_y] == .x$value[!.x$is_y]))
) %>%
unlist()
if (is.null(K)) {
K <- min(
length(unique(ppc_data(y, yrep)$rep_id)) + 1,
length(pit)
)
K <- nrow(yrep) + 1
}
} else {
inform("'pit' specified so ignoring 'y', and 'yrep' if specified.")
Expand All @@ -632,20 +632,22 @@ ppc_pit_ecdf <- function(y,
ggplot() +
aes(
x = 1:K / K,
y = ecdf(pit)(seq(0, 1, length.out = K)) - (plot_diff == TRUE) * 1:K / K,
y = ecdf(pit)(seq(0, 1, length.out = K)) -
(plot_diff == TRUE) * seq(0, 1, length.out = K),
color = "y"
) +
geom_step(show.legend = FALSE) +
geom_step(aes(
y = lims$upper[-1] / N - (plot_diff == TRUE) * 1:K / K,
y = lims$upper[-1] / N - (plot_diff == TRUE) * seq(0, 1, length.out = K),
color = "yrep"
), show.legend = FALSE) +
),
linetype = 2, show.legend = FALSE) +
geom_step(aes(
y = lims$lower[-1] / N - (plot_diff == TRUE) * 1:K / K,
y = lims$lower[-1] / N - (plot_diff == TRUE) * seq(0, 1, length.out = K),
color = "yrep"
), show.legend = FALSE) +
yaxis_title(FALSE) +
xaxis_title(FALSE) +
),
linetype = 2, show.legend = FALSE) +
labs(y = "ECDF", x = "PIT") +
yaxis_ticks(FALSE) +
scale_color_ppc() +
bayesplot_theme_get()
Expand All @@ -671,10 +673,13 @@ ppc_pit_ecdf_grouped <-
if (is.null(pit)) {
pit <- ppc_data(y, yrep, group) %>%
group_by(.data$y_id) %>%
dplyr::group_map(~ mean(.x$value[.x$is_y] >= .x$value[!.x$is_y])) %>%
dplyr::group_map(
~ mean(.x$value[.x$is_y] > .x$value[!.x$is_y]) +
runif(1, max = mean(.x$value[.x$is_y] == .x$value[!.x$is_y]))
) %>%
unlist()
if (is.null(K)) {
K <- length(unique(ppc_data(y, yrep)$rep_id)) + 1
K <- nrow(yrep) + 1
}
} else {
inform("'pit' specified so ignoring 'y' and 'yrep' if specified.")
Expand Down Expand Up @@ -723,13 +728,14 @@ ppc_pit_ecdf_grouped <-
geom_step(aes(
y = .data$lims_upper - (plot_diff == TRUE) * .data$x,
color = "yrep"
), show.legend = FALSE) +
),
linetype = 2, show.legend = FALSE) +
geom_step(aes(
y = .data$lims_lower - (plot_diff == TRUE) * .data$x,
color = "yrep"
), show.legend = FALSE) +
xaxis_title(FALSE) +
yaxis_title(FALSE) +
),
linetype = 2, show.legend = FALSE) +
labs(y = "ECDF", x = "PIT") +
yaxis_ticks(FALSE) +
bayesplot_theme_get() +
facet_wrap("group") +
Expand Down
Binary file modified R/sysdata.rda
Binary file not shown.

0 comments on commit 7e00829

Please sign in to comment.