R/sher_monte_carlo.R

In jlovegren0/aermod: Work with AERMOD POST files

#' Generate 99.9 percentile value of total exceedances (10,000 Monte Carlo trials)
#' for State Health Effects review, given table of exceedances.
#' 
#' Generates various statistics used in TCEQ Toxicology Tier III RFC reviews.
#' If the impact matrix covers more than one year of met data, only the first year is used.
#'
#' @name sher_monte_carlo
#' @param lambda A lambda object generated by the \code{sher_lambda} function
#' @param nhr A proposed limit on hr/yr of operatoin.
#' @return A table of 99.9% percentile exceedance figures for different window sizes and ESL thresholds (2xESL, 4xESL)
#' @export
sher_monte_carlo <- function(lambda,nhr,nsim=10000L){
	exc_tbl <- lambda$exc_tbl
	SMPL <- function( x , blocklen = 24L ){ 
		starting_idxs <- 1 : length(x)
		x <- c( x , x[1:( blocklen - 1)] )
		p_tbl <- purrr::map_int( starting_idxs , function(.x) sum( x[ .x : ( .x + blocklen - 1 ) ] ) ) %>% 
			tibble::tibble( nexc = . ) %>% dplyr::group_by(nexc) 
		res <- p_tbl %>% dplyr::tally()
		names(res)[2] <- paste0('n_',blocklen)
		res
	}
	ptbl <- function( x = "e2" )
	{
		probs <- tibble::tibble( nexc = 0:24 , dummy = 0 ) %>%
		dplyr::left_join( SMPL( dplyr::pull(exc_tbl,x) ,blocklen = 1L) ) %>% 
		dplyr::left_join( SMPL( dplyr::pull(exc_tbl,x) ,blocklen = 2L) ) %>% 
		dplyr::left_join( SMPL( dplyr::pull(exc_tbl,x) ,blocklen = 3L) ) %>% 
		dplyr::left_join( SMPL( dplyr::pull(exc_tbl,x) ,blocklen = 4L) ) %>% 
		dplyr::left_join( SMPL( dplyr::pull(exc_tbl,x) ,blocklen = 6L) ) %>% 
		dplyr::left_join( SMPL( dplyr::pull(exc_tbl,x) ,blocklen = 12L) ) %>% 
		dplyr::left_join( SMPL( dplyr::pull(exc_tbl,x) ,blocklen = 24L) ) %>% dplyr::select(-dummy) %>% suppressMessages
		probs[ is.na(probs) ] <- 0L
		probs %>% dplyr::mutate( lval = x ) %>% dplyr::select( 1 , 9 , 2:8 )
	}
	pbz <- dplyr::bind_rows( ptbl("e1"), ptbl( "e2" ) , ptbl( "e4" ), ptbl("e10") )
	nxc <- function(lx,blocksz){
		((0:24) * stats::rmultinom( nsim , size = floor( nhr / blocksz ), dplyr::pull(dplyr::filter(pbz,lval==lx),paste0('n_',blocksz)))) %>%
			colSums %>% quantile( 0.999 ) %>% as.integer
	}
	res_frame <- tibble::tibble(dummy=1,lx=c('e1','e2','e4','e10')) %>% 
		dplyr::left_join(tibble::tibble(dummy=1,blocksz=c(1,2,3,4,6,12,24))) %>%
		dplyr::select(-dummy) %>% dplyr::mutate( nexc = 0 ) %>% suppressMessages
	res_frame$nexc <- purrr::map_int( 1:nrow(res_frame) , ~ nxc( res_frame$lx[.] , res_frame$blocksz[.] ) )
	res_frame %<>% tidyr::pivot_wider( names_from = lx , values_from = nexc )
	worst <- res_frame %>% dplyr::summarize( L1 = max(e1) , L2 = max(e2) , L4 = max(e4) , L10 = max(e10) )
	list( nhr = nhr , results = res_frame , worst_lambda = worst)
}