mrc-ide / drjacoby

Flexible Markov chain monte carlo via reparameterization
https://mrc-ide.github.io/drjacoby/
Other
12 stars 5 forks source link

Parameter summary #103

Open bobverity opened 2 years ago

bobverity commented 2 years ago

Function to return common summaries of all parameters from the output, e.g. ML, MAP, posterior median, credible intervals etc. Think about if/how these should also be part of the standard output.

pwinskill commented 2 years ago

Some thing like:

par_summary <- function(mcmc, lower_quantile = 0.025, upper_quantile = 0.975){
  samples <- mcmc$output %>% 
    dplyr::filter(.data$phase == "sampling")

  mode <- samples %>%
    slice_max(order_by = (.data$logprior + .data$loglikelihood), n = 1, with_ties = FALSE) %>%
    select(-c(.data$chain, .data$phase, .data$iteration, .data$logprior, .data$loglikelihood)) %>%
    tidyr::pivot_longer(-c(), names_to = "Parameter", values_to = "mode")

 par_summary <- samples %>%
    dplyr::select(-c(.data$chain, .data$phase, .data$iteration, .data$logprior, .data$loglikelihood)) %>%
    tidyr::pivot_longer(cols = -c(), names_to = "Parameter", values_to = "Value") %>%
    dplyr::group_by(.data$Parameter) %>%
    dplyr::summarise(
      mean = mean(.data$Value),
      median = median(.data$Value),
      lower_CrI = quantile(.data$Value, lower_quantile),
      upper_CrI = quantile(.data$Value, upper_quantile),
    ) %>%
    dplyr::left_join(mode, by = "Parameter") %>%
    dplyr::select(.data$Parameter, .data$mean, .data$mode, .data$median, .data$lower_CrI, .data$upper_CrI)

 return(par_summary)
}

This works on the drjacoby_output object without any down-sampling at present