tidy_conditional_event_prior_sampler.R

library(tidyverse)

### nuber of reruns ###
nsims <- 10000


## parameter table defining the parameter space --------------------
params <- tibble(param_name = c("NNIG", "NPUE", "NANC", "TDIV", "TMIG",  "MIGRPROBN2P", "MIGRPROBP2N"),
                 log10min   = c(rep(3,3), 1, 1, rep(-5, 2)),
                 log10max   = c(rep(log10(2.5e5), 3), rep(4, 2), rep(-2, 2)),
                 units = c(rep("haploid inds", 3), rep("generations", 2), rep("m", 2)),
                 log10mid = (log10min+log10max)/2)

## the actual sampler ----------------------------------------------
param_table <- tibble(NNIG = round(10^runif(nsims, min = params$log10min[1], max = params$log10max[1])),
                      NPUE = round(10^runif(nsims, min = params$log10min[2], max = params$log10max[2])),
                      NANC = round(10^runif(nsims, min = params$log10min[3], max = params$log10max[3])),
                      TDIV = runif(nsims, min = params$log10min[4], max = params$log10max[4]),
                      #TMIG = round(TDIV* runif(nsims, min = 0, max = 1)), # uniform [0-1] * max
                      TMIG = runif(nsims, min = params$log10min[5], max = TDIV),                        # the coditional part
                      MIGRPROBN2P = 10^runif(nsims, min = params$log10min[6], max = params$log10max[6]),
                      MIGRPROBP2N = 10^runif(nsims, min = params$log10min[7], max = params$log10max[7])) %>%
  mutate(TDIV = round(10^TDIV),
         TMIG = round(10^TMIG))


## summary table for the vertical lines ------------------------------------------
sumtab <- param_table %>%
  gather(key = 'param_name', value = 'val') %>% 
  group_by(param_name) %>%
  summarise(mean = mean(val),
            median = median(val))


# distribution among the parameter space of the two conditional parameters -------
param_table %>%
  ggplot() +
  geom_point(aes(x = TDIV, y = TMIG))+
  # geom_hex(aes(x = TDIV, y = TMIG),bins=80)+
  scale_x_log10()+
  scale_y_log10()+
  scale_fill_viridis_c()

# the old historgams -------------------------------------------------------------
param_table %>%
  gather(key = 'param_name', value = 'val') %>% 
  ggplot(aes(x=val)) +
  #geom_density() +
  geom_histogram(color='black',fill='lightgray')+
  geom_vline(data = sumtab, aes(xintercept = mean), color = "red",size = .5)+
  geom_vline(data = sumtab, aes(xintercept = median), color = "red",size = .5,linetype=2)+
  geom_vline(data = params,aes(xintercept = 10^log10mid), color = "black",size = .5)+
  scale_x_log10()+
  facet_wrap(param_name~.,scales = 'free')