R/tv_upper_bound.R

In pierrejacob/montecarlodsm: Gibbs sampler for Dempster-Shafer inference in Categorical distributions

# dempsterpolytope - Gibbs sampler for Dempster's inference in Categorical distributions 
# 
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.

#'@rdname tv_upper_bound
#'@title Turns meeting times into upper bounds on TV to stationarity
#'@description
#' Given a vector of meeting times, e.g. generated by \code{\link{sample_meeting_times}}, a choice of lag L
#' and an iteration t, returns
#' mean(pmax(0,ceiling((meetingtimes-L-t)/L))).
#'@param meetingtimes vector of meeting times, such as generated by \code{\link{sample_meeting_times}}.
#'@param lag lag between the chains.
#'@param t iteration at which to compute the TV upper bound.
#'@return A non-negative value representing an 
#'estimated upper bound between the chain at iteration t and the limiting distribution
#'of the chain.
#'@examples
#'\dontrun{
#'nrep <- 100
#'meeting_times <- sapply(1:nrep, function(irep) sample_meeting_times(counts = c(3,2,0,1)))
#'tmax <- floor(max(meeting_times)*1.2)
#'ubounds <- sapply(1:tmax, function(t) tv_upper_bound(meeting_times, 1, t))
#'plot(x = 1:tmax, y = ubounds, type = 'l', xlab = "iteration", ylab = "TV upper bounds")
#'}
#'@export
tv_upper_bound <- function(meetingtimes, L, t){
  return(mean(pmax(0,ceiling((meetingtimes-L-t)/L))))
}