diff --git a/BayesianTools/NAMESPACE b/BayesianTools/NAMESPACE index 6ca782d..71021bb 100644 --- a/BayesianTools/NAMESPACE +++ b/BayesianTools/NAMESPACE @@ -68,6 +68,7 @@ export(getPossibleSamplerTypes) export(getPredictiveDistribution) export(getPredictiveIntervals) export(getSample) +export(getSampleDistance) export(getVolume) export(likelihoodAR1) export(likelihoodIidNormal) diff --git a/Test/Distance.R b/BayesianTools/R/Distance.R similarity index 77% rename from Test/Distance.R rename to BayesianTools/R/Distance.R index 4f06e8d..e36b747 100644 --- a/Test/Distance.R +++ b/BayesianTools/R/Distance.R @@ -11,6 +11,9 @@ #' BH = Bhattacharyya distance #' D = normalized Euclidean distance between mean and standard deviation of sample and target. Not symmetric. Target is sample 1. This was described in eq. 10 in Laloy, E., and J. A. Vrugt. 2012. High-dimensional posterior exploration of hydrologic models using multiple-try DREAM(ZS) and high-performance computing. Water Resour. Res. 48(1) #' +#' @example /inst/examples/getSampleDistanceHelp.R +#' @export +#' getSampleDistance <- function(sample1, sample2, type = "KL"){ if(type == "KL"){ @@ -18,7 +21,7 @@ getSampleDistance <- function(sample1, sample2, type = "KL"){ x = FNN::KL.dist(sample1, sample2, k=10) out = mean(x) # FH: no idea if the mean is a good idea. KL.dist returns one value per cluster size, I don't know which value is best chosen - } else if(type == "BH" ){ + } else if(type == "BH" | type == "BHs" ){ Sigma1 = cov(sample1) Sigma2 = cov(sample2) mu1 = colMeans(sample1) @@ -27,10 +30,14 @@ getSampleDistance <- function(sample1, sample2, type = "KL"){ # The following code is copied from package fpc aggregatesigma <- (Sigma1+Sigma2)/2 - d1 <- mahalanobis(mu1,mu2,aggregatesigma)/8 - d2 <- log(det(as.matrix(aggregatesigma))/sqrt(det(as.matrix(Sigma1))* - det(as.matrix(Sigma2))))/2 + #d1 <- mahalanobis(mu1,mu2,aggregatesigma,tol=1e-20)/8 # tol (tolerance) argument gets passed to solve(). Avoids crashes if there are very small values in the covariance matrices + d1 <- mahalanobis(mu1,mu2,aggregatesigma,tol=.Machine$double.xmin)/8 # tol (tolerance) argument gets passed to solve(). Avoids crashes if there are very small values in the covariance matrices + d2 <- log(det(as.matrix(aggregatesigma))/sqrt(det(as.matrix(Sigma1))*det(as.matrix(Sigma2))))/2 + if(type == "BH"){ out <- d1+d2 + } else if (type == "BHs"){ + out <- c(d1, d2) + } # end fpc @@ -42,7 +49,7 @@ getSampleDistance <- function(sample1, sample2, type = "KL"){ sd2 = apply(sample2, 2, sd) mu1 = colMeans(sample1) mu2 = colMeans(sample2) - dev = sum( ((mu1 - mu2)/sd1)^2 + ((mu1 - mu2)/sd1)^2 ) + dev = sum( ((mu1 - mu2)/sd1)^2 + ((sd1 - sd2)/sd1)^2 ) out = sqrt( 1/(2*length(mu1)) * dev ) @@ -72,23 +79,7 @@ getSampleDistance <- function(sample1, sample2, type = "KL"){ #' -library(mvtnorm) - -sigma <- matrix(c(4,2,2,3), ncol=2) -X <- rmvnorm(n=2500, mean=c(1,2), sigma=sigma) -X = rbind(X,X) -Y <- rmvnorm(n=2500, mean=c(5,2), sigma=sigma) -Y = rbind(Y,X) - -FNN::KL.dist(X, Y, k=10) - -X = getSample(out, start = 1000, end = 2000) -Y = getSample(out, start = 2000, end = 3000) - -getSampleDistance(X,Y, type = "KL") -getSampleDistance(X,Y, type = "BH") -getSampleDistance(X,Y, type = "D") diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index d740eb5..4947125 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -1,113 +1,821 @@ #' SMC sampler -#' @author Florian Hartig -#' @description Sequential Monte Carlo Sampler -#' @param bayesianSetup either an object of class bayesianSetup created by \code{\link{createBayesianSetup}} (recommended), or a log target function -#' @param initialParticles initial particles - either a draw from the prior, provided as a matrix with the single parameters as columns and each row being one particle (parameter vector), or a numeric value with the number of desired particles. In this case, the sampling option must be provided in the prior of the BayesianSetup. +#' +#' @description A Sequential Monte Carlo (SMC) Sampler with differential evolution updating to avoid particle depletion +#' +#' @param bayesianSetup either an object of class bayesianSetup created by \code{\link{createBayesianSetup}} (recommended), or a log target function +#' @param initialParticles initial particles - either a draw from the prior, provided as a matrix with the single parameters as columns and each row being one particle (parameter vector), or a numeric value with the number of desired particles. In this case, the sampling option must be provided in the prior of the BayesianSetup. #' @param iterations number of iterations -#' @param resampling logical, specifies whether new particles should be created at each iteration -#' @param resamplingSteps how many resampling (MCMC) steps between the iterations -#' @param proposal optional, proposal class -#' @param adaptive logical, should the covariance of the proposal be adapted during sampling? -#' @param proposalScale scaling factor for the proposal generation. Can be adapted if there is too much / too little rejection +#' @param exponents series of exponents to build the intermediate distributions +#' @param lastMutateSteps how many resampling (MCMC) steps after the SMC iterations to increase sample diversity. +#' @param proposal optional proposal class +#' @param x Parameter to generate the exponential sequence for building intermediary distributions. Default value from Jeremiah et al. (2012) +#' @param m Parameter to generate the exponential sequence for building intermediary distributions. Default value from Jeremiah et al. (2012) +#' @param sampling Which algorithm to use for particle (re)sampling. Options are "multinomial" (default), "residual" and "systematic" +#' @param ess.limit Threshold value (as a fraction between 0 and 1) of the effective sample compared to the number of particles at which resampling is done. To resample at each step, use a value > 1. +#' @param lastResample Iteration (starting from the end) at which particle resampling is forced. To deactivate this, set to a value < 0 +#' +#' @param resampling deprecated, use resamplingSteps to modify resampling +#' @param resamplingSteps how many resampling (MCMC) steps between the SMC iterations. A value of 20 or higher is recommended. +#' @param mutate.method MCMC method for the resampling. Either "Metropolis" (default) or "DE" +#' @param adaptive should the covariance of the proposal be adapted during re-sampling +#' @param proposalScale scaling factor for the proposal generation. Can be adapted if there is too much / too little rejection. Experience from MCMC sampling suggests an optimal value of 2.38/sqrt(2*d) (where d is the number of parameters), and first experiments suggest that this value is also appropriate for SMC sampling. +#' @param pars.lower Optional: vector containing the minimum values of calibration parameters. Used to initialise the proposal function +#' @param pars.upper Optional: vector containing the maximum values of calibration parameters. Used to initialise the proposal function +#' @param b parameter for snooker update in DE, set by default to 10e-4. As b implies a certain scale for the parameters, b should be either scaled to the parameter uncertainty, or parameters should be scaled. +#' @param diagnostics an optional function with diagnostics that are calculated on the particles during each SMC iteration +#' #' @details The sampler can be used for rejection sampling as well as for sequential Monte Carlo. For the former case set the iterations to one. +#' @note The SMC currently assumes that the initial particle is sampled from the prior. If a better initial estimate of the posterior distribution is available, this the sampler should be modified to include this. Currently, however, this is not included in the code, so the appropriate adjustments have to be done by hand. +#' +#' @author Florian Hartig, Matthias Speich #' -#' @note The SMC currently assumes that the initial particle is sampled from the prior. If a better initial estimate of the posterior distribution is available, this the sampler should be modified to include this. Currently, however, this is not included in the code, so the appropriate adjustments have to be done by hand. #' @export #' @example /inst/examples/SMCHelp.R -smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, resampling = T, resamplingSteps = 2, proposal = NULL, adaptive = T, proposalScale = 0.5){ - - if(resamplingSteps < 1) stop("SMC error, resamplingSteps can't be < 1") - - setup <- checkBayesianSetup(bayesianSetup) - +smcSampler <- function(bayesianSetup, + initialParticles = 1000, + iterations = 10, + resampling = NULL, + resamplingSteps = 2, + lastMutateSteps = 5, + proposal = NULL, + exponents = NULL, + adaptive = T, + proposalScale = 0.5, + x=3.11, + m=7E-08, + sampling="multinomial", + ess.limit=0.5, + ess.factor = 0.95, + lastResample = 1, + pars.lower=NULL, + pars.upper=NULL, + mutate.method ="Metropolis", + b=1e-04, + diagnostics = NULL, + reference=NULL){ + + # SETUP STEPS + + if(!is.null(resampling)) message("resampling is deprecated, use resamplingSteps to control resampling") + + # Timing: normally done in mcmcRun.R. Timing is included in this function for the purpose of testing the SMC function. + # Starting the clock + ptm <- proc.time() + + #if(resamplingSteps < 1) stop("SMC error, resamplingSteps can't be < 1") + info = list() - info$resamplingAcceptance = matrix(nrow = iterations, ncol = resamplingSteps) - info$survivingParticles = rep(NA, iterations) - - - if(inherits(initialParticles, "numeric")){ + + # The number of iterations is not known at the beginning. Therefore, output vectors (with one value for each iteration) + # are made large enough that they will probably not need to be grown (growing arrays can be slow). + if(resamplingSteps > 0) info$resamplingAcceptance = as.data.frame(matrix(nrow = 10000, ncol = resamplingSteps)) #Using DF instead of matrix because number of iterations may change due to adaptive algorithm + info$survivingParticles = rep(NA, 10000) + info$ess.vec <- info$exponents <- vector("numeric", 10000) + info$diagnostics <- diag.end <- list() + lastAccept <- vector("numeric", lastMutateSteps) + if(ess.limit < 0) stop("ess.limit can't be < 0") + + setup <- checkBayesianSetup(bayesianSetup) + + ### InitialParticles + #TODO documentation + + if(class(initialParticles)[1] == "numeric"){ initialParticles = bayesianSetup$prior$sampler(initialParticles) + importanceDensity = bayesianSetup$prior$density } - + if(class(initialParticles) == "matrix"){ + importanceDensity = bayesianSetup$prior$density + } + if(class(initialParticles) == "list"){ + importanceDensity = initialParticles$density + if(class(initialParticles$particles) == "numeric"){ + initialParticles = as.matrix(initialParticles$particles,ncol=1) + } else if(class(initialParticles$particles) == "data.frame"){ + initialParticles = as.matrix(initialParticles$particles) + } else if(class(initialParticles$particles) == "matrix"){ + initialParticles <- initialParticles$particles + } + } + if (any(is.infinite(setup$prior$density(initialParticles)))) stop("initialParticles outside prior range") - + particles <- initialParticles - rejectionRate = 0 - particleSize = nrow(initialParticles) - - acceptanceTarget = round(particleSize / 2) - + ess.limit.abs = round(particleSize) * ess.limit + + weights <- oldweights <- oldInter <- rep(0, particleSize) + posterior = matrix(nrow = particleSize, ncol = 3) - numPar <- ncol(initialParticles) - - if (is.null(proposal)) proposalGenerator = createProposalGenerator(rep(40,numPar)) - + + if(mutate.method %in% c("Metropolis", "adaptive")){ + if (is.null(proposal)) proposalGenerator = createProposalGenerator(rep(40,numPar)) + } + usedUp = 0 - - for (i in 1:iterations){ - - posterior = setup$posterior$density(particles, returnAll = T) - - likelihoodValues <- posterior[,2] - - # idea - adjust (1/iterations) such that always approx 30% of particles are maintain - #level = sort(likelihoodValues)[acceptanceTarget] - #best = likelihoodValues - - ll = likelihoodValues - max(likelihoodValues, na.rm = T) - - llCutoff = sort(ll)[acceptanceTarget] - - - - relativeL = exp(likelihoodValues - max(likelihoodValues, na.rm = T))^(1/iterations) - - sel = sample.int(n=length(likelihoodValues), size = length(likelihoodValues), replace = T, prob = relativeL) - - info$survivingParticles[i] = length(unique(sel)) - - particles = particles[sel,] - + + #### TEST + #proposalScale <- runif(nrow(particles), min = 1e-06, max = 1) + Z <- particles + + + # WEIGHT SETUP + + # Define an exponential sequence of beta parameters, i.e. the exponents to weight the likelihood against prior. + # Idea and parameter values from Jeremiah et al. (2012), Environ Modell Softw + # if(is.null(exponents)){ + # iters <- seq(0,iterations) + # exponents <- m * (iters*200/iterations)^x + # exponents <- pmin(1,exponents) + # exponents <- exponents[2:length(exponents)] + # } + + # Initial value for exponent: here the first value of the series is taken. + # If the situation ESS < E* never occurs, sampler execution is the same as + # with the non-adaptive algorithm + estar <- round(0.2*particleSize) + oldExp <- 0 + #curExp <- exponents[1] + curExp <- 0 # For fully adaptive algorithm: initial value of curExp does not matter, but must be < 1 + + icount <- 1 # Iterations counter + + # All particles are given equal weight at initialization + weights[1:length(weights)] <- log(1/particleSize) + oldweights <- weights + + # Initial importance distribution + importanceValues <- importanceDensity(particles) + #print("Calculating initial posteriors") + posteriorValues <- setup$posterior$density(particles) + #print("Initial posteriors done") + #oldInter <- importanceValues + oldExp <- 0 + + if(!is.null(pars.lower) & !is.null(pars.lower)){ + sds <- 0.1 * (pars.upper - pars.lower) + } else{ + sds <- rep(40,numPar) + } + if (is.null(proposal)){ proposalGenerator = createProposalGenerator(sds)} + + lastIteration <- FALSE + + ## Estimate minimum and maximum value of each parameter (from initial population) + min.pars <- apply(particles,2,min) + max.pars <- apply(particles,2,max) + + + #for (i in 1:iterations){ + # while(icount <= length(exponents)){ + #while(curExp <= 1 & !lastIteration){ # For fully adaptive algorithm + # if(curExp == 1){lastIteration <- TRUE} + while(curExp < 1){ + + ### TEMPORARY DEBUG + if(sum(is.nan(posteriorValues)) > 0 | sum(is.na(posteriorValues)) > 0 | sum(is.infinite(posteriorValues)) > 0){ + print("NaN / NA / Infinite in importanceValues") + if(sum(is.nan(posteriorValues)) > 0){ + print(c("NaN indices", which(is.nan(posteriorValues)))) + print("Particles with NaN:") + print(particles[is.nan(posteriorValues),]) + } + + if(sum(is.na(posteriorValues)) > 0){ + print(c("NA indices", which(is.na(posteriorValues)))) + print("Particles with NA:") + print(particles[is.na(posteriorValues),]) + } + + if(sum(is.infinite(posteriorValues)) > 0){ + print(c("Infinite indices", which(is.infinite(posteriorValues)))) + print("Particles with Infinite:") + print(particles[is.infinite(posteriorValues),]) + } + } + + # Update particle min/max (if a particle is more "extreme" than previously recorded) + cur.min <- apply(particles,2,min) + cur.max <- apply(particles,2,max) + min.pars <- pmin(min.pars, cur.min) + max.pars <- pmax(max.pars, cur.max) + if (numPar == 1) particles = matrix(particles, ncol = 1) - - if (resampling == T){ - - if (adaptive == T){ - proposalGenerator = updateProposalGenerator(proposalGenerator, particles) + + # Using a while loop instead of for because in the adaptive algorithm, the number of iterations may increase + + # Mutate + + # 1) When to do the resampling + # 2) HOW - currently adaptive Metropolis, could also do DEzs step + + # Terminology issue: in literature, "resampling" generally describes the replication of particles with high likelihood, i.e. the previous + # step. This step here is generally referred to as "mutate" or "move". + # It may be necessary to rename the parameters "resampling" and "resamplingSteps", as well as info$resamplingAcceptance + + + # mutate.out <- mutate(setup = setup, particles = particles, proposalGenerator = proposalGenerator, posteriorValues = posteriorValues, importanceDensity = importanceDensity, method = mutate.method, steps = resamplingSteps, proposalScale = proposalScale, adaptive = adaptive, b=b) + # particles <- mutate.out$particles + # posteriorValues <- mutate.out$posteriorValues + # importanceValues <- mutate.out$importanceValues + # info$resamplingAcceptance[(icount),] <- mutate.out$acceptance + + # Reweighting + + # idea - adjust (1/iterations) such that always approx 30% of particles are maintain + # level = sort(likelihoodValues)[acceptanceTarget] + # best = likelihoodValues + + # Jeremia is using beta sequence + + # https://github.com/florianhartig/BayesianTools/issues/23 + + + ## Update weights + # Create intermediary distribution + ess <- 1 / sum(exp(2 * weights)) + oldInter <- oldExp * posteriorValues + (1-oldExp) * importanceValues + + # TODO 8.3.26 - what's up with this alternative call using ess.factor, which was also present ? For the moment, I erased this option from the code. + # ess.factor was in the main function call + #inter.out <- beta.search(ess=ess, target.ess = (ess * ess.factor), posteriorValues = posteriorValues, importanceValues = importanceValues, oldInter = oldInter, curWeights = weights, curExp = curExp) + inter.out <- beta.search(ess=ess, + target.ess = (ess.limit.abs-1), # TODO 8.3.26 why -1? + posteriorValues = posteriorValues, + importanceValues = importanceValues, + oldInter = oldInter, + curWeights = weights, + curExp = curExp) + curExp <- inter.out$newExp + weights <- inter.out$weights + interDist <- inter.out$interDist + ess <- inter.out$ess + doResample <- inter.out$doResample + + + info$ess.vec[icount] <- ess + info$exponents[icount] <- curExp + if(abs(1-curExp) < 1e-03){ + curExp <- 1 + } + + # Resampling step + + # Determine if resampling is necessary - if yes, resample with probability given by the weight + # Resample also on the last iteration, or at the iteration (starting from the end) given by the parameter lastResample (output is based + # on the location of particles in parameter space, the weights are not considered in the output) + + if(ess < ess.limit.abs | icount == (length(exponents) - lastResample) | doResample){ + + oldExp <- curExp + + sel = resample(weights, method = sampling) + + particles = particles[sel,] + + #print(apply(particles, 2, sd)) + #print(apply(particles, 2, sd)/apply(particles, 2, mean)) + + posteriorValues <- posteriorValues[sel] + importanceValues <- importanceValues[sel] + info$survivingParticles[icount] = length(unique(sel)) + + # Set all weights equal + weights[1:length(weights)] <- log(1/particleSize) + # Normalize (log-)weights so that the sum (of non-logs) equals 1 + weights <- weights - BayesianTools:::logSumExp(weights) + + ## Check diversity of particles, and require a mutate step if necessary + # Scale particles + #print(c("nrow unique particles", nrow(unique(particles)))) + #print(c("curExp", curExp)) + #if(nrow(unique(particles)) < 1000){ + scaled.particles <- t(apply(particles, 1, scale.particles, min.pars = min.pars, max.pars = max.pars)) + #diversity <- rao.div(scaled.particles = scaled.particles) + #print(c("Diversity", diversity)) + #} + + if(resamplingSteps > 0){ + mutate.out <- mutate(setup = setup, particles = particles, proposalGenerator = proposalGenerator, posteriorValues = posteriorValues, importanceDensity = importanceDensity, method = mutate.method, steps = resamplingSteps, proposalScale = proposalScale, adaptive = adaptive, b=b, min.pars = min.pars, max.pars = max.pars, Z = NULL) + particles <- mutate.out$particles + posteriorValues <- mutate.out$posteriorValues + importanceValues <- mutate.out$importanceValues + info$resamplingAcceptance[(icount),] <- mutate.out$acceptance + if(length(proposalScale) > 1){ + proposalScale <- mutate.out$proposalScale } - - for (j in 1:resamplingSteps){ - particlesProposals = proposalGenerator$returnProposalMatrix(particles, scale = proposalScale) - - jumpProb <- exp(setup$posterior$density(particlesProposals) - likelihoodValues[sel])^(i/iterations) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) - - accepted <- jumpProb > runif(length(jumpProb), 0 ,1) - - rejectionRate = rejectionRate + sum(accepted) - - particles[accepted, ] = particlesProposals[accepted, ] - - + + #plot(density(particles[,1]), xlim=c(0.314, 0.681), main = c(curExp, nrow(unique(particles)))) + + #Z <- rbind(Z, particles) + + # while(nrow(unique(particles)) < nrow(particles) * 0.5){ + # mutate.out <- mutate(setup = setup, particles = particles, proposalGenerator = proposalGenerator, posteriorValues = posteriorValues, importanceDensity = importanceDensity, method = mutate.method, steps = resamplingSteps, proposalScale = proposalScale, adaptive = adaptive, b=b, min.pars = min.pars, max.pars = max.pars, Z = NULL) + # + # particles <- mutate.out$particles + # posteriorValues <- mutate.out$posteriorValues + # importanceValues <- mutate.out$importanceValues + # info$resamplingAcceptance[(icount),] <- mutate.out$acceptance + # if(length(proposalScale) > 1){ + # proposalScale <- mutate.out$proposalScale + # } + # + # plot(density(particles[,1]), xlim=c(0.314, 0.681), main = c("oi", curExp, nrow(unique(particles)))) + # } + } + } + + ## Mutate step if required because of nuerical issues or low particle diversity + # if(doResample){ + # mutate.out <- mutate(setup = setup, particles = particles, proposalGenerator = proposalGenerator, posteriorValues = posteriorValues, importanceDensity = importanceDensity, method = mutate.method, steps = resamplingSteps, proposalScale = proposalScale, adaptive = adaptive, b=b) + # particles <- mutate.out$particles + # posteriorValues <- mutate.out$posteriorValues + # importanceValues <- mutate.out$importanceValues + # info$resamplingAcceptance[(icount),] <- mutate.out$acceptance + # } + + oldweights <- weights + info$ess.vec[icount] <- ess + + icount <- icount + 1 + oldExp <- curExp + + if(!is.null(diagnostics)) info$diagnostics[[icount-1]] <- diagnostics(particles, reference) + } + + if (numPar == 1) particles = matrix(particles, ncol = 1) + + # Last resampling step, so that particles are distributed according to target distribution + sel = resample(weights, method = sampling) + particles = particles[sel,] + info$survivingParticles[(icount-1)] <- length(unique(sel)) + + #print("last") + #print(apply(particles, 2, sd)/apply(particles, 2, mean)) + + if(resamplingSteps > 0){ + for(mutateStep in 1:lastMutateSteps){ + # Last mutation, to increase diversity between particles + mutate.out <- mutate(setup = setup, particles = particles, proposalGenerator = proposalGenerator, posteriorValues = posteriorValues, importanceDensity = importanceDensity, method = mutate.method, steps = 1, proposalScale = proposalScale, adaptive = adaptive, b=b, min.pars = min.pars, max.pars = max.pars, Z = NULL) + particles <- mutate.out$particles + posteriorValues <- mutate.out$posteriorValues + importanceValues <- mutate.out$importanceValues + lastAccept[mutateStep] <- mutate.out$acceptance + if(!is.null(diagnostics)) diag.end[[mutateStep]] <- diagnostics(particles, reference) + + #plot(density(particles[,1]), xlim=c(0.314, 0.681)) + + #Z <- rbind(Z, particles) + } + } - info$rejectionRate = rejectionRate / (iterations * resamplingSteps) + # Timing: normally done in mcmcRun.R. Timing is included in this function for the purpose of testing the SMC function. + # Stopping the clock + elapsed.time <- proc.time() - ptm + + # Trim info objects (as their size = number of iterations were not known beforehand, they + # were allocated a large vector/matrix). Trim up to (icount-1), as icount got incremented + # at the end of the last iteration. + if(resamplingSteps > 0) info$resamplingAcceptance <- info$resamplingAcceptance[1:(icount-1),] + info$survivingParticles <- info$survivingParticles[1:(icount-1)] + info$exponents <- info$exponents[1:(icount-1)] + #info$diagnostics <- c(info$diagnostics[[1:(icount-1)]], diag.end) + info$diagnostics <- c(info$diagnostics, diag.end) + info$ess.vec <- info$ess.vec[1:(icount-1)] + if(resamplingSteps > 0) info$lastAccept <- lastAccept + info$elapsed.time <- elapsed.time + + settings = list(initialParticles = initialParticles, proposalGenerator = proposalGenerator) + out = list( setup = setup, - initialParticles = initialParticles, + settings = settings, particles = particles, - posteriorValues = posterior, - proposalGenerator = proposalGenerator, + posterior = posterior, info = info ) - + class(out) <- c("smcSampler", "bayesianOutput") return(out) - + } + + +#' Residual Resampling +#' +#' @keywords internal +residualResampling <- function(weights){ + weights <- exp(weights) + # Define number of replications for each particle (integer), as well as residuals + nrep <- weights * length(weights) + nrep.int <- floor(nrep) + residuals <- nrep - nrep.int + + # Resample remaining particles based on residuals + length.missing <- length(weights) - sum(nrep.int) + missing <- vector("numeric", length.missing) + missing <- sample.int(n=length(weights), size=length.missing, replace = TRUE, prob = residuals) + + # Return the indices of particles to be resampled + new.parts <- rep(1:length(weights), nrep.int) + new.parts <- c(new.parts, missing) + return(new.parts) +} + +systematicResampling <- function(weights){ + weights <- exp(weights) + # Reorder weights in increasing order + rank.weights <- rank(weights, ties.method = "first") + sort.weights <- sort(weights) + cumu.weights <- c(0,cumsum(sort.weights)[1:(length(weights)-1)]) + n.parts <- length(weights) # Number of particles + + u <- runif(n=1, min=0, max=1/n.parts) + U <- vector("numeric", n.parts) + new.parts <- vector("numeric", n.parts) + new.parts.sort <- vector("numeric", n.parts) + for(i in 1:n.parts){ + U[i] <- ((i-1)/n.parts) + u + + + new.parts.sort[i] <- tail(which(cumu.weights <= U[i]), 1) + # Match sorted ranks back to original order of weights + new.parts[i] <- which(rank.weights==new.parts.sort[i]) + } + + return(new.parts) +} + + +#' Auxiliary function for adaptive SMC algorithm +#' +#' +#' @keywords internal +resample <- function(weights, method = "multinomial"){ + + particleSize = length(weights) + + if(method == "multinomial"){ + sel <- sample.int(n=particleSize, size = particleSize, replace = T, prob = exp(weights)) + } else if(method == "residual"){ + sel <- residualResampling(weights) + } else if(method == "systematic"){ + sel <- systematicResampling(weights) + } else{ + stop("Invalid string for resampling argument") + } + return(sel) +} + + +#' Auxiliary function for adaptive SMC algorithm +#' +#' +#' @keywords internal +beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldInter, curWeights, curExp, tol=1){ + # A function to dynamically set the next exponent to build the next intermediary distribution. + # Uses the bisection method. Following Jasra et al. (2011), Scand J Statist, doi: 10.1111/j.1467-9469.2010.00723.x + + # Initial exponent - set to 1 (maximum possible value) + tryDiff <- 1 - curExp + a <- curExp + b <- 1 + try.ess <- target.ess + (100*tol) # Dummy initial value that is sure to be different from target value + tryExp <- 0 + + + while(abs(try.ess-target.ess) > tol & tryExp < 1 & abs(b-a) > 1e-10){ + tryExp <- (a+b) * 0.5 + + + tryDist <- tryExp * posteriorValues + (1-tryExp) * importanceValues + + tryWeights <- curWeights + (tryDist - oldInter) + # Normalize (log-)weights so that the sum (of non-logs) equals 1 + + ### DEBUG + if(sum(is.na(tryWeights)) > 0){ + print("DEBUG") + print("tryWeights") + print(tryWeights) + print("tryDist") + print(tryDist) + print("oldInter") + print(oldInter) + print(c("curExp", curExp)) + print(c("tryExp", tryExp)) + print(c("a", a)) + print(c("b", b)) + print("posteriorValues") + print(posteriorValues) + print("importanceValues") + print(importanceValues) + + } + + tryWeights <- tryWeights - BayesianTools:::logSumExp(tryWeights) + try.ess <- 1 / sum(exp(2 * tryWeights)) + + if(try.ess - target.ess > tol & !any(is.infinite(tryWeights))){ + # Greater ESS than desired -> choose a *larger* exponent in the next iteration (efficiency) + a <- tryExp + } else if(try.ess - target.ess < -tol | any(is.infinite(tryWeights))){ + # Smaller ESS than desired -> choose a *smaller* exponent in the next iteration (stability) + # Also includes a failsafe in the case of infinite weights (numerical issue) + b <- tryExp + } + } + + + + if(any(is.infinite(tryWeights)) | is.infinite(try.ess)){ + doResample <- TRUE + tryExp <- curExp + tryWeights <- curWeights + tryDist <- oldInter + try.ess <- ess + } else{ + doResample <- FALSE + } + + out <- list(newExp=tryExp, weights=tryWeights, interDist=tryDist, ess=try.ess, doResample = doResample) + return(out) +} + + +#' Auxiliary function for particle mutation +#' +#' +#' @keywords internal +mutate <- function(setup, particles, proposalGenerator, posteriorValues, importanceDensity, method, steps, proposalScale, adaptive = TRUE, b=1E-04, min.pars=NULL, max.pars=NULL, Z = NULL){ + + if(is.vector(particles)){particles = matrix(particles, ncol = 1)} + acceptance <- vector("numeric", length=steps) + importanceValues <- importanceDensity(particles) + + + if(length(proposalScale) == 1){ + scale.factors <- rep(proposalScale, nrow(particles)) + } else { + scale.factors <- proposalScale + scaled.particles <- t(apply(particles, 1, scale.particles, min.pars = min.pars, max.pars = max.pars)) + } + + if(is.null(Z)){ + Z <- particles + } + + r_extra <- vector("numeric", length = nrow(particles)) + + + if(adaptive){ + proposalGenerator = updateProposalGenerator(proposalGenerator, particles) + } + + for(j in 1:steps){ + + if(method=="Metropolis"){ + particlesProposals = proposalGenerator$returnProposalMatrix(particles, scale = proposalScale) + } else if(method=="DE"){ + + particlesProposals <- particles + + for(part in 1:nrow(particles)){ + + if(runif(1) < 0.1) { + ## Snooker update + particlesOld <- particles + particleDiff12 <- rep(0,ncol(particles)) + particleDiff23 <- rep(0,ncol(particles)) + particleDiff13 <- rep(0,ncol(particles)) + newParts <- rep(0L,3) + + # Sample 3 other particles + #while(all(particleDiff12==0) | all(particleDiff23==0) | all(particleDiff13==0) | part %in% newParts){ + while(all(particleDiff12==0) | all(particleDiff23==0) | all(particleDiff13==0)){ + # The sampled particles might be identical (especially after resampling). Therefore, it is checked whether + # the difference between particles is non-zero for at least one parameter. If the particles are identical, + # the sampled particles are discarded and two new particles are sampled. + # Also, the sampled particles should not include the current particle. + newParts <- sample.int(n=nrow(particles),size=3, replace=FALSE) + particleDiff12 <- Z[newParts[1],] - Z[newParts[2],] + particleDiff23 <- Z[newParts[2],] - Z[newParts[3],] + particleDiff13 <- Z[newParts[1],] - Z[newParts[3],] + } + + x_z <- particles[part,] - Z[newParts[3],] + D2 <- max(sum(x_z*x_z), 1.0e-300) + projdiff <- sum((Z[newParts[1],] -Z[newParts[2],]) * x_z)/D2 # inner_product of difference with x_z / squared norm x_z + particlesProposals[part,] <- particles[part,] + runif(1, min=1.2, max=2.2) * projdiff * x_z + + x_z <- particlesProposals[part,] - Z[newParts[3],] + D2prop <- max(sum(x_z*x_z), 1.0e-300) + r_extra[part] <- ((ncol(particles)-1)/2) * (log(D2prop) - log(D2)) + + #if(part %in% 1:10){ + # print(c("Part snooker", part, r_extra[part])) + #} + + } else{ + ## No snooker + particlesOld <- particles + particleDiff <- rep(0,ncol(particles)) + newParts <- rep(0L,2) + + # Sample 2 other particles + #while(all(particleDiff==0) | part %in% newParts){ + while(all(particleDiff==0)){ + # The sampled particles might be identical (especially after resampling). Therefore, it is checked whether + # the difference between particles is non-zero for at least one parameter. If the particles are identical, + # the sampled particles are discarded and two new particles are sampled. + # Also, the sampled particles should not include the current particle. + newParts <- sample.int(n=nrow(particles),size=2, replace=FALSE) + particleDiff <- Z[newParts[1],] - Z[newParts[2],] + } + numPar <- ncol(particles) + randVector <- runif(numPar,-b,b) + particlesProposals[part,] <- particles[part,] + (particleDiff * scale.factors[part]) + randVector + r_extra[part] <- 0 + + # if(part %in% 1:6){ + # print(c("Part ", part)) + # print("Particle") + # print(particlesOld[part,]) + # print("Proposal") + # print(particlesProposals[part,]) + # } + #if(part %in% 1:10){ + # print(c("Part NO snooker", part, r_extra[part])) + #} + } + } + } + + + proposalPosteriors <- setup$posterior$density(particlesProposals) + proposalImportance <- importanceDensity(particlesProposals) + + #print(c("Priors", head(setup$prior$density(particlesProposals), 10))) + #print(c("Proposals post", head(proposalPosteriors, 10))) + #print(c("Particles post", head(posteriorValues, 10))) + + #print(c("Infinite", sum(is.infinite(proposalPosteriors)))) + #print(c("Infinite prior", sum(is.infinite(setup$prior$density(particlesProposals))))) + + #print(c("Minus Infinite", sum(proposalPosteriors== -Inf))) + #print(c("Minus Infinite prior", sum(setup$prior$density(particlesProposals)== -Inf))) + + ### TODO: why is jumpProb not simply the ratio of posteriors, i.e. why include the prior? + # This has no influence on the current experiments (with flat priors), but needs to be clarified before other uses. + + #jumpProb <- exp(proposalPosteriors - posteriorValues) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) + jumpProb <- exp(proposalPosteriors - posteriorValues + r_extra) + + accepted <- jumpProb > runif(length(jumpProb), 0 ,1) + + if(sum(is.na(accepted)) > 1){ + print("sum is.na accepted", sum(is.na(accepted))) + print("proposalPosteriors") + print(proposalPosteriors) + print("posteriorValues") + print(posteriorValues) + print("r_extra") + print(r_extra) + + } + + particles[accepted, ] = particlesProposals[accepted, ] + posteriorValues[accepted] <- proposalPosteriors[accepted] + importanceValues[accepted] <- proposalImportance[accepted] + acceptance[j] <- sum(accepted)/nrow(particles) + + ## Adaptation of scaling factors following Fearnhead & Taylor (2013) + if(length(proposalScale) > 1){ + + ### Calculate ESJD for each proposal + ##Scale particles + scaled.proposals <- t(apply(particlesProposals, 1, scale.particles, min.pars = min.pars, max.pars = max.pars)) + + + #print("scaled.particles") + #print(head(scaled.particles)) + + #print("scaled.proposals") + #print(head(scaled.proposals)) + + + #esjd <- vector("numeric", nrow(particles)) + esjd <- vector("numeric", sum(r_extra == 0)) + loop.ind <- 0 + ind.nosnooker <- which(r_extra == 0) + for(part in 1:nrow(particles)){ + if(r_extra[part] == 0) { + loop.ind <- loop.ind + 1 + esjd[loop.ind] <- dist(rbind(scaled.particles[part,], scaled.proposals[part,])) * pmin(1,jumpProb[part]) + } + #print(rbind(scaled.particles[part,], scaled.proposals[part,])) + #if(r_extra[part] != 0){print("Snooker")} + #print(c("dist, jumpProb, post, prop",dist(rbind(scaled.particles[part,], scaled.proposals[part,])), jumpProb[part], posteriorValues[part], proposalPosteriors[part])) + } + + #print(c("length esjd", length(esjd))) + + esjd[is.infinite(esjd)] <- 0 + #print(c("esjd", head(esjd))) + #print(c("max esjd", max(esjd))) + + # Normalize + esjd <- esjd/sum(esjd) + #print(c("esjd norm", head(esjd))) + #print(c("max esjd norm", max(esjd))) + + # Resample and add noise + sel <- resample(log(esjd), "systematic") + #print(c("sel", head(sel))) + #print(c("unique sel", unique(sel))) + if(length(unique(sel)) <= 5){ + scale.factors <- scale.factors #+ rnorm(length(scale.factors), 0, 0.015) + } else{ + scale.factors <- scale.factors[ind.nosnooker[sel]] #+ rnorm(length(scale.factors), 0, 0.015) + } + + # Draw scale factors from normal distribution + mean.sf <- mean(scale.factors) + sd.sf <- sd(scale.factors) + + scale.factors <- rnorm(n = nrow(particles), mean = mean.sf, sd = max(0.05,sd.sf)) + + scale.factors[scale.factors <= 0] <- 1E-06 + + } + + #print(head(scale.factors, 100)) + #plot(density(scale.factors), xlim=c(0,1)) + + } + #} + if(length(proposalScale) == 1){ + out <- list(particles=particles, posteriorValues = posteriorValues, importanceValues = importanceValues, acceptance = acceptance) + } else{ + out <- list(particles=particles, posteriorValues = posteriorValues, importanceValues = importanceValues, acceptance = acceptance, proposalScale = scale.factors) + } + + return(out) +} + + +#' Scale particles +#' +#' @keywords internal +scale.particles <- function(particle.row, min.pars, max.pars){ + scaled.particles <- (particle.row - min.pars) / (max.pars - min.pars) + return(scaled.particles) +} + + +#' Calculate rao diversity +#' +#' @keywords internal +rao.div <- function(scaled.particles){ + # Determine unique particles, their sum and pairwise difference + unique.particles <- unique(scaled.particles) + unique.count <- vector("numeric", nrow(unique.particles)) + for(i in 1:nrow(unique.particles)){ + unique.count[i] <- sum(apply(scaled.particles, 1, function(x, cur.unique){ + if(sum(x-cur.unique)==0){ + return(TRUE) + } else { + return(FALSE) + } + }, cur.unique = unique.particles[i,])) + } + + unique.count <- unique.count/nrow(scaled.particles) + + total.dist <- 0 + + for(i in 1:(nrow(unique.particles)-1)){ + for(j in (i+1):nrow(unique.particles)){ + #total.dist <- total.dist + dist(unique.particles[c(i,j),]) * unique.count[i] * unique.count[j] + total.dist <- total.dist + sum(abs(unique.particles[i,] - unique.particles[j,])) * unique.count[i] * unique.count[j] + } + } + + return(total.dist) +} + + + +# Auxiliary function to calculate Metropolis ratio when log-likelihood values are very large + +# largeRatio <- function(a,b){ +# # Determine the larger of the two numbers +# +# } diff --git a/BayesianTools/R/classMcmcSampler.R b/BayesianTools/R/classMcmcSampler.R index a00c51e..e05d517 100644 --- a/BayesianTools/R/classMcmcSampler.R +++ b/BayesianTools/R/classMcmcSampler.R @@ -140,7 +140,8 @@ summary.mcmcSampler <- function(object, printCorrelation = "auto", ...){ psf <- FALSE mcmcsampler <- sampler$settings$sampler - runtime <- sampler$settings$runtime[3] + runtime <- sampler$info$runtime[3] + correlations <- round(cor(getSample(sampler)),3) chain <- getSample(sampler, parametersOnly = T, coda = T, ...) # chain <- getSample(sampler, parametersOnly = T, coda = T) diff --git a/BayesianTools/R/classMcmcSamplerList.R b/BayesianTools/R/classMcmcSamplerList.R index e3a9512..2d5f6ee 100644 --- a/BayesianTools/R/classMcmcSamplerList.R +++ b/BayesianTools/R/classMcmcSamplerList.R @@ -38,7 +38,7 @@ summary.mcmcSamplerList <- function(object, ...){ mcmcsampler <- sampler[[1]]$settings$sampler runtime <- 0 - for(i in 1:length(sampler)) runtime <- runtime+sampler[[i]]$settings$runtime[3] + for(i in 1:length(sampler)) runtime <- runtime+sampler[[i]]$info$runtime[3] correlations <- round(cor(getSample(sampler)),3) diff --git a/BayesianTools/R/mcmcRun.R b/BayesianTools/R/mcmcRun.R index 2e4f7d1..27066d6 100644 --- a/BayesianTools/R/mcmcRun.R +++ b/BayesianTools/R/mcmcRun.R @@ -248,9 +248,11 @@ runMCMC <- function(bayesianSetup , sampler = "DEzs", settings = NULL){ mcmcSampler$settings = settings } - mcmcSampler$settings$runtime = mcmcSampler$settings$runtime + proc.time() - ptm + mcmcSampler$info$sessionInfo = utils::sessionInfo() + + mcmcSampler$info$runtime = mcmcSampler$info$runtime + proc.time() - ptm if(is.null(settings$message) || settings$message == TRUE){ - message("runMCMC terminated after ", mcmcSampler$settings$runtime[3], "seconds") + message("runMCMC terminated after ", mcmcSampler$info$runtime[3], "seconds") } return(mcmcSampler) } @@ -310,64 +312,84 @@ applySettingsDefault<-function(settings=NULL, sampler = "DEzs", check = FALSE){ if(!settings$sampler %in% getPossibleSamplerTypes()$BTname) stop("trying to set values for a sampler that does not exist") - - mcmcDefaults <- list(startValue = NULL, - iterations = 10000, - burnin = 0, - thin = 1, - consoleUpdates = 100, - parallel = NULL, - message = TRUE) - - #### Metropolis #### - if(settings$sampler %in% c("AM", "DR", "DRAM", "Metropolis")){ - - defaultSettings <- c(mcmcDefaults, list(optimize = T, - proposalGenerator = NULL, - adapt = F, - adaptationInterval = 500, - adaptationNotBefore = 3000, - DRlevels = 1 , - proposalScaling = NULL, - adaptationDepth = NULL, - temperingFunction = NULL, - proposalGenerator = NULL, - gibbsProbabilities = NULL)) - - if (settings$sampler %in% c("AM", "DRAM")) defaultSettings$adapt <- TRUE - if (settings$sampler %in% c("DR", "DRAM")) defaultSettings$DRlevels <- 2 + if (settings$sampler == "AM") { + defaultSettings <- applySettingsDefault(sampler = "Metropolis") + defaultSettings$adapt <- TRUE } - #### DE Family #### - if(settings$sampler %in% c("DE", "DEzs")){ - defaultSettings <- c(mcmcDefaults, list(eps = 0, - currentChain = 1, - blockUpdate = list("none", - k = NULL, - h = NULL, - pSel = NULL, - pGroup = NULL, - groupStart = 1000, - groupIntervall = 1000) - )) - - if (settings$sampler == "DE"){ - defaultSettings$f <- -2.38 # TODO CHECK - - } - - if (settings$sampler == "DEzs"){ - defaultSettings$f <- 2.38 - defaultSettings <- c(defaultSettings, list(Z = NULL, - zUpdateFrequency = 1, - pSnooker = 0.1, - pGamma1 = 0.1, - eps.mult =0.2, - eps.add = 0)) - } - + if (settings$sampler == "DR") { + defaultSettings <- applySettingsDefault(sampler = "Metropolis") + defaultSettings$DRlevels <- 2 + } + + if (settings$sampler == "DRAM") { + defaultSettings <- applySettingsDefault(sampler = "Metropolis") + defaultSettings$adapt <- TRUE + defaultSettings$DRlevels <- 2 } + if (settings$sampler == "Metropolis"){ + defaultSettings = list(startValue = NULL, + iterations = 10000, + optimize = T, + proposalGenerator = NULL, + consoleUpdates=100, + burnin = 0, + thin = 1, + parallel = NULL, + adapt = T, + adaptationInterval= 500, + adaptationNotBefore = 3000, + DRlevels = 1 , + proposalScaling = NULL, + adaptationDepth = NULL, + temperingFunction = NULL, + proposalGenerator = NULL, + gibbsProbabilities = NULL, + currentChain = 1, + message = TRUE) + } + + if (settings$sampler == "DE"){ + defaultSettings = list(startValue = NULL, + iterations = 10000, + burnin = 0, + thin = 1, + eps = 0, + consoleUpdates = 100, + currentChain = 1, + parallel = F, + f = -2.38, # TODO CHECK + burnin = 0, + eps = 0, + consoleUpdates = 100, + blockUpdate = list("none", k = NULL, h = NULL, pSel = NULL, pGroup = NULL, + groupStart = 1000, groupIntervall = 1000), + message = TRUE) + } + + if (settings$sampler == "DEzs"){ + defaultSettings = list(startValue = NULL, + iterations = 10000, + Z = NULL, + pSnooker = 0.1, + burnin = 0, + thin = 1, + f = 2.38, # TODO CHECK + eps = 0, + pGamma1 = 0.1, + eps.mult =0.2, + eps.add = 0, + consoleUpdates = 100, + currentChain = 1, + parallel = NULL, + zUpdateFrequency = 1, + blockUpdate = list("none", k = NULL, h = NULL, pSel = NULL, pGroup = NULL, + groupStart = 1000, groupIntervall = 1000), + message = TRUE) + } + + #### DREAM Family #### if(settings$sampler %in% c("DREAM", "DREAMzs")){ @@ -410,18 +432,30 @@ applySettingsDefault<-function(settings=NULL, sampler = "DEzs", check = FALSE){ #### SMC #### if (settings$sampler == "SMC"){ - defaultSettings = list(iterations = 10, - resampling = T, - resamplingSteps = 2, - proposal = NULL, - adaptive = T, - proposalScale = 0.5, - initialParticles = 1000 + defaultSettings = list( initialParticles = 1000, + iterations = 10, + resampling = T, + resamplingSteps = 2, + lastMutateSteps = 5, + proposal = NULL, + exponents = NULL, + adaptive = T, + proposalScale = 0.5, + x=3.11, + m=7E-08, + sampling="multinomial", + ess.limit=NULL, + ess.factor = 0.95, + lastResample = 1, + pars.lower=NULL, + pars.upper=NULL, + mutate.method ="Metropolis", + b=1e-04, + diagnostics = NULL, + reference=NULL ) } - - ## CHECK DEFAULTS if(check){ @@ -439,9 +473,7 @@ applySettingsDefault<-function(settings=NULL, sampler = "DEzs", check = FALSE){ } defaultSettings$nrChains = 1 - defaultSettings$runtime = 0 - defaultSettings$sessionInfo = utils::sessionInfo() - + nam = names(defaultSettings) for (i in 1:length(defaultSettings)){ diff --git a/BayesianTools/inst/examples/SMCHelp.R b/BayesianTools/inst/examples/SMCHelp.R index 8270740..5704dc3 100644 --- a/BayesianTools/inst/examples/SMCHelp.R +++ b/BayesianTools/inst/examples/SMCHelp.R @@ -1,94 +1,19 @@ -## Example for the use of SMC -# First we need a bayesianSetup - SMC makes most sense if we can for demonstration, -# we'll write a function that puts out the number of model calls - -MultiNomialNoCor <- generateTestDensityMultiNormal(sigma = "no correlation") - -parallelLL <- function(parMatrix){ - print(paste("Calling likelihood with", nrow(parMatrix), "parameter combinations")) - out = apply(parMatrix, 1, MultiNomialNoCor) - return(out) -} - -bayesianSetup <- createBayesianSetup(likelihood = parallelLL, lower = rep(-10, 3), +# creating a simple test BayesianSetup +bayesianSetup <- createBayesianSetup(likelihood = generateTestDensityMultiNormal(), + lower = rep(-10, 3), upper = rep(10, 3), parallel = "external") -# Defining settings for the sampler -# First we use the sampler for rejection sampling -settings <- list(initialParticles = 1000, iterations = 1, resampling = FALSE) - -# Running the sampler -out1 <- runMCMC(bayesianSetup = bayesianSetup, sampler = "SMC", settings = settings) -#plot(out1) - - -# Now for sequential Monte Carlo -settings <- list(initialParticles = 100, iterations = 5, resamplingSteps = 1) -out2 <- runMCMC(bayesianSetup = bayesianSetup, sampler = "SMC", settings = settings) -#plot(out2) - -\dontrun{ - -## Example for starting a new SMC run with results from a previous SMC run - -# Generate example data (time series) -# x1 and x2 are predictory, yObs is the response -t <- seq(1, 365) -x1 <- (sin( 1 / 160 * 2 * pi * t) + pi) * 5 -x2 <- cos( 1 / 182.5 * 1.25 * pi * t) * 12 - -# the model -mod <- function(par, t1 = 1, tn = 365) { - par[1] * x1[t1:tn] + par[2] * x2[t1:tn] -} - -# the true parameters -par1 <- 1.65 -par2 <- 0.75 -yObs <- mod(c(par1, par2)) + rnorm(length(x1), 0, 2) - -# split the time series in half -plot(yObs ~ t) -abline(v = 182, col = "red", lty = 2) - -# First half of the data -ll_1 <- function(x, sum = TRUE) { - out <- dnorm(mod(x, 1, 182) - yObs[1:182], 0, 2, log = TRUE) - if (sum == TRUE) sum(out) else out -} - -# Fit the first half of the time series -# (e.g. fit the model to the data soon as you collect the data) -setup_1 <- createBayesianSetup(ll_1, lower = c(-10, -10), upper = c(10, 10)) -settings_1 <- list(initialParticles = 1000) -out_1 <- runMCMC(setup_1, "SMC", settings_1) -summary(out_1) - - -# Second half of the data -ll_2 <- function(x, sum = TRUE) { - out <- dnorm(mod(x, 183, 365) - yObs[183:365], 0, 2, log = TRUE) - if (sum == TRUE) sum(out) else out -} - -# Fit the second half of the time series -# (e.g. fit the model to the data soon as you collect the data) -setup_2 <- createBayesianSetup(ll_2, lower = c(-10, -10), upper = c(10, 10)) - -# This is the important step, we use the final particles from the -# previous SMC run to initialize the new SMC run -settings_2 <- list(initialParticles = out_1$particles) -out_2 <- runMCMC(setup_2, "SMC", settings_2) -summary(out_2) - -par_pred <- apply(out_2$particles, 2, median) -pred <- mod(par_pred) -plotTimeSeries(yObs, pred) - -} - - +# to directly call the smcSampler, use +SMCout <- smcSampler(bayesianSetup, initialParticles = 100, iterations = 3) +# plotting results +marginalPlot(SMCout, type = "d", singlePanel = T) +# for options, see +?smcSampler +# as for all MCMCs, you can also run the SMC via the general runMCMC function +SMCout<-runMCMC(bayesianSetup, sampler = "SMC", + settings = list(nrChains =2, initialParticles = 100, iterations = 3)) +plot(SMCout) diff --git a/BayesianTools/inst/examples/getSampleDistanceHelp.R b/BayesianTools/inst/examples/getSampleDistanceHelp.R new file mode 100644 index 0000000..6f06340 --- /dev/null +++ b/BayesianTools/inst/examples/getSampleDistanceHelp.R @@ -0,0 +1,37 @@ +library(mvtnorm) +library(BayesianTools) + +sigma <- matrix(c(4,2,2,3), ncol=2) +X <- rmvnorm(n=2500, mean=c(1,2), sigma=sigma) +X = rbind(X,X) +Y <- rmvnorm(n=2500, mean=c(5,2), sigma=sigma) +Y = rbind(Y,X) + +FNN::KL.dist(X, Y, k=10) + +getSampleDistance(X,Y, type = "KL") +getSampleDistance(X,Y, type = "BH") +getSampleDistance(X,Y, type = "D") + +# the distance functions are independent of the scale of the parameters, see e.g. + +scale = 2 +getSampleDistance(data.frame(rnorm(1000), rnorm(1000, sd = scale)), + data.frame(rnorm(1000), rnorm(1000, mean = scale, sd = scale)), type = "BH") +getSampleDistance(data.frame(rnorm(1000), rnorm(1000, sd = scale)), + data.frame(rnorm(1000), rnorm(1000, mean = scale, sd = scale)), type = "D") + + +# To get an idea about how the distance functions approach zero for a perfect sampler +# use the following code where we sample 2 distributions directly from their +# rnorm function + +getD <- function(n){ + x1 = matrix(rnorm(n), ncol = 10) + x2 = matrix(rnorm(n), ncol = 10) + getSampleDistance(x1, x2, type = "D" ) +} + +test = seq(1000, 30000, 1000) +plot(test, sapply(test, getD), ylim = c(0,0.4)) +abline(h=0, col = "red") diff --git a/BayesianTools/man/beta.search.Rd b/BayesianTools/man/beta.search.Rd new file mode 100644 index 0000000..9811291 --- /dev/null +++ b/BayesianTools/man/beta.search.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/SMC.R +\name{beta.search} +\alias{beta.search} +\title{Auxiliary function for adaptive SMC algorithm} +\usage{ +beta.search( + ess, + target.ess, + posteriorValues, + importanceValues, + oldInter, + curWeights, + curExp, + tol = 1 +) +} +\description{ +Auxiliary function for adaptive SMC algorithm +} +\keyword{internal} diff --git a/BayesianTools/man/getSampleDistance.Rd b/BayesianTools/man/getSampleDistance.Rd new file mode 100644 index 0000000..2dc346d --- /dev/null +++ b/BayesianTools/man/getSampleDistance.Rd @@ -0,0 +1,64 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Distance.R +\name{getSampleDistance} +\alias{getSampleDistance} +\title{calculate distance between samples} +\usage{ +getSampleDistance(sample1, sample2, type = "KL") +} +\arguments{ +\item{sample1}{a matrix. If the distance function is not symmetric, this is assumed to be a sample from the target (reference) distribution} + +\item{sample2}{a matrix. If the distance function is not symmetric, this is assumed to be a sample from the "other" distribution} + +\item{type}{the distance function} +} +\description{ +calculate distance between samples +} +\details{ +Currently, the following distance functions are implemented + +KL = Kullback-Leibler Divergence +BH = Bhattacharyya distance +D = normalized Euclidean distance between mean and standard deviation of sample and target. Not symmetric. Target is sample 1. This was described in eq. 10 in Laloy, E., and J. A. Vrugt. 2012. High-dimensional posterior exploration of hydrologic models using multiple-try DREAM(ZS) and high-performance computing. Water Resour. Res. 48(1) +} +\examples{ +library(mvtnorm) +library(BayesianTools) + +sigma <- matrix(c(4,2,2,3), ncol=2) +X <- rmvnorm(n=2500, mean=c(1,2), sigma=sigma) +X = rbind(X,X) +Y <- rmvnorm(n=2500, mean=c(5,2), sigma=sigma) +Y = rbind(Y,X) + +FNN::KL.dist(X, Y, k=10) + +getSampleDistance(X,Y, type = "KL") +getSampleDistance(X,Y, type = "BH") +getSampleDistance(X,Y, type = "D") + +# the distance functions are independent of the scale of the parameters, see e.g. + +scale = 2 +getSampleDistance(data.frame(rnorm(1000), rnorm(1000, sd = scale)), + data.frame(rnorm(1000), rnorm(1000, mean = scale, sd = scale)), type = "BH") +getSampleDistance(data.frame(rnorm(1000), rnorm(1000, sd = scale)), + data.frame(rnorm(1000), rnorm(1000, mean = scale, sd = scale)), type = "D") + + +# To get an idea about how the distance functions approach zero for a perfect sampler +# use the following code where we sample 2 distributions directly from their +# rnorm function + +getD <- function(n){ + x1 = matrix(rnorm(n), ncol = 10) + x2 = matrix(rnorm(n), ncol = 10) + getSampleDistance(x1, x2, type = "D" ) +} + +test = seq(1000, 30000, 1000) +plot(test, sapply(test, getD), ylim = c(0,0.4)) +abline(h=0, col = "red") +} diff --git a/BayesianTools/man/mutate.Rd b/BayesianTools/man/mutate.Rd new file mode 100644 index 0000000..9056bf6 --- /dev/null +++ b/BayesianTools/man/mutate.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/SMC.R +\name{mutate} +\alias{mutate} +\title{Auxiliary function for particle mutation} +\usage{ +mutate( + setup, + particles, + proposalGenerator, + posteriorValues, + importanceDensity, + method, + steps, + proposalScale, + adaptive = TRUE, + b = 1e-04, + min.pars = NULL, + max.pars = NULL, + Z = NULL +) +} +\description{ +Auxiliary function for particle mutation +} +\keyword{internal} diff --git a/BayesianTools/man/rao.div.Rd b/BayesianTools/man/rao.div.Rd new file mode 100644 index 0000000..60800f8 --- /dev/null +++ b/BayesianTools/man/rao.div.Rd @@ -0,0 +1,12 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/SMC.R +\name{rao.div} +\alias{rao.div} +\title{Calculate rao diversity} +\usage{ +rao.div(scaled.particles) +} +\description{ +Calculate rao diversity +} +\keyword{internal} diff --git a/BayesianTools/man/resample.Rd b/BayesianTools/man/resample.Rd new file mode 100644 index 0000000..5409123 --- /dev/null +++ b/BayesianTools/man/resample.Rd @@ -0,0 +1,12 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/SMC.R +\name{resample} +\alias{resample} +\title{Auxiliary function for adaptive SMC algorithm} +\usage{ +resample(weights, method = "multinomial") +} +\description{ +Auxiliary function for adaptive SMC algorithm +} +\keyword{internal} diff --git a/BayesianTools/man/residualResampling.Rd b/BayesianTools/man/residualResampling.Rd new file mode 100644 index 0000000..55bec9e --- /dev/null +++ b/BayesianTools/man/residualResampling.Rd @@ -0,0 +1,12 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/SMC.R +\name{residualResampling} +\alias{residualResampling} +\title{Residual Resampling} +\usage{ +residualResampling(weights) +} +\description{ +Residual Resampling +} +\keyword{internal} diff --git a/BayesianTools/man/scale.particles.Rd b/BayesianTools/man/scale.particles.Rd new file mode 100644 index 0000000..e37389f --- /dev/null +++ b/BayesianTools/man/scale.particles.Rd @@ -0,0 +1,12 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/SMC.R +\name{scale.particles} +\alias{scale.particles} +\title{Scale particles} +\usage{ +\method{scale}{particles}(particle.row, min.pars, max.pars) +} +\description{ +Scale particles +} +\keyword{internal} diff --git a/BayesianTools/man/smcSampler.Rd b/BayesianTools/man/smcSampler.Rd index a1da4bb..e927013 100644 --- a/BayesianTools/man/smcSampler.Rd +++ b/BayesianTools/man/smcSampler.Rd @@ -8,11 +8,25 @@ smcSampler( bayesianSetup, initialParticles = 1000, iterations = 10, - resampling = T, + resampling = NULL, resamplingSteps = 2, + lastMutateSteps = 5, proposal = NULL, + exponents = NULL, adaptive = T, - proposalScale = 0.5 + proposalScale = 0.5, + x = 3.11, + m = 7e-08, + sampling = "multinomial", + ess.limit = 0.5, + ess.factor = 0.95, + lastResample = 1, + pars.lower = NULL, + pars.upper = NULL, + mutate.method = "Metropolis", + b = 1e-04, + diagnostics = NULL, + reference = NULL ) } \arguments{ @@ -22,121 +36,70 @@ smcSampler( \item{iterations}{number of iterations} -\item{resampling}{logical, specifies whether new particles should be created at each iteration} +\item{resampling}{deprecated, use resamplingSteps to modify resampling} -\item{resamplingSteps}{how many resampling (MCMC) steps between the iterations} +\item{resamplingSteps}{how many resampling (MCMC) steps between the SMC iterations. A value of 20 or higher is recommended.} -\item{proposal}{optional, proposal class} +\item{lastMutateSteps}{how many resampling (MCMC) steps after the SMC iterations to increase sample diversity.} -\item{adaptive}{logical, should the covariance of the proposal be adapted during sampling?} +\item{proposal}{optional proposal class} -\item{proposalScale}{scaling factor for the proposal generation. Can be adapted if there is too much / too little rejection} -} -\description{ -Sequential Monte Carlo Sampler -} -\details{ -The sampler can be used for rejection sampling as well as for sequential Monte Carlo. For the former case set the iterations to one. -} -\note{ -The SMC currently assumes that the initial particle is sampled from the prior. If a better initial estimate of the posterior distribution is available, this the sampler should be modified to include this. Currently, however, this is not included in the code, so the appropriate adjustments have to be done by hand. -} -\examples{ -## Example for the use of SMC -# First we need a bayesianSetup - SMC makes most sense if we can for demonstration, -# we'll write a function that puts out the number of model calls +\item{exponents}{series of exponents to build the intermediate distributions} -MultiNomialNoCor <- generateTestDensityMultiNormal(sigma = "no correlation") +\item{adaptive}{should the covariance of the proposal be adapted during re-sampling} -parallelLL <- function(parMatrix){ - print(paste("Calling likelihood with", nrow(parMatrix), "parameter combinations")) - out = apply(parMatrix, 1, MultiNomialNoCor) - return(out) -} +\item{proposalScale}{scaling factor for the proposal generation. Can be adapted if there is too much / too little rejection. Experience from MCMC sampling suggests an optimal value of 2.38/sqrt(2*d) (where d is the number of parameters), and first experiments suggest that this value is also appropriate for SMC sampling.} -bayesianSetup <- createBayesianSetup(likelihood = parallelLL, lower = rep(-10, 3), - upper = rep(10, 3), parallel = "external") +\item{x}{Parameter to generate the exponential sequence for building intermediary distributions. Default value from Jeremiah et al. (2012)} -# Defining settings for the sampler -# First we use the sampler for rejection sampling -settings <- list(initialParticles = 1000, iterations = 1, resampling = FALSE) - -# Running the sampler -out1 <- runMCMC(bayesianSetup = bayesianSetup, sampler = "SMC", settings = settings) -#plot(out1) +\item{m}{Parameter to generate the exponential sequence for building intermediary distributions. Default value from Jeremiah et al. (2012)} +\item{sampling}{Which algorithm to use for particle (re)sampling. Options are "multinomial" (default), "residual" and "systematic"} -# Now for sequential Monte Carlo -settings <- list(initialParticles = 100, iterations = 5, resamplingSteps = 1) -out2 <- runMCMC(bayesianSetup = bayesianSetup, sampler = "SMC", settings = settings) -#plot(out2) +\item{ess.limit}{Threshold value (as a fraction between 0 and 1) of the effective sample compared to the number of particles at which resampling is done. To resample at each step, use a value > 1.} -\dontrun{ +\item{lastResample}{Iteration (starting from the end) at which particle resampling is forced. To deactivate this, set to a value < 0} -## Example for starting a new SMC run with results from a previous SMC run +\item{pars.lower}{Optional: vector containing the minimum values of calibration parameters. Used to initialise the proposal function} -# Generate example data (time series) -# x1 and x2 are predictory, yObs is the response -t <- seq(1, 365) -x1 <- (sin( 1 / 160 * 2 * pi * t) + pi) * 5 -x2 <- cos( 1 / 182.5 * 1.25 * pi * t) * 12 +\item{pars.upper}{Optional: vector containing the maximum values of calibration parameters. Used to initialise the proposal function} -# the model -mod <- function(par, t1 = 1, tn = 365) { - par[1] * x1[t1:tn] + par[2] * x2[t1:tn] -} - -# the true parameters -par1 <- 1.65 -par2 <- 0.75 -yObs <- mod(c(par1, par2)) + rnorm(length(x1), 0, 2) +\item{mutate.method}{MCMC method for the resampling. Either "Metropolis" (default) or "DE"} -# split the time series in half -plot(yObs ~ t) -abline(v = 182, col = "red", lty = 2) +\item{b}{parameter for snooker update in DE, set by default to 10e-4. As b implies a certain scale for the parameters, b should be either scaled to the parameter uncertainty, or parameters should be scaled.} -# First half of the data -ll_1 <- function(x, sum = TRUE) { - out <- dnorm(mod(x, 1, 182) - yObs[1:182], 0, 2, log = TRUE) - if (sum == TRUE) sum(out) else out +\item{diagnostics}{an optional function with diagnostics that are calculated on the particles during each SMC iteration} } - -# Fit the first half of the time series -# (e.g. fit the model to the data soon as you collect the data) -setup_1 <- createBayesianSetup(ll_1, lower = c(-10, -10), upper = c(10, 10)) -settings_1 <- list(initialParticles = 1000) -out_1 <- runMCMC(setup_1, "SMC", settings_1) -summary(out_1) - - -# Second half of the data -ll_2 <- function(x, sum = TRUE) { - out <- dnorm(mod(x, 183, 365) - yObs[183:365], 0, 2, log = TRUE) - if (sum == TRUE) sum(out) else out +\description{ +A Sequential Monte Carlo (SMC) Sampler with differential evolution updating to avoid particle depletion } - -# Fit the second half of the time series -# (e.g. fit the model to the data soon as you collect the data) -setup_2 <- createBayesianSetup(ll_2, lower = c(-10, -10), upper = c(10, 10)) - -# This is the important step, we use the final particles from the -# previous SMC run to initialize the new SMC run -settings_2 <- list(initialParticles = out_1$particles) -out_2 <- runMCMC(setup_2, "SMC", settings_2) -summary(out_2) - -par_pred <- apply(out_2$particles, 2, median) -pred <- mod(par_pred) -plotTimeSeries(yObs, pred) - +\details{ +The sampler can be used for rejection sampling as well as for sequential Monte Carlo. For the former case set the iterations to one. } +\note{ +The SMC currently assumes that the initial particle is sampled from the prior. If a better initial estimate of the posterior distribution is available, this the sampler should be modified to include this. Currently, however, this is not included in the code, so the appropriate adjustments have to be done by hand. +} +\examples{ +# creating a simple test BayesianSetup +bayesianSetup <- createBayesianSetup(likelihood = generateTestDensityMultiNormal(), + lower = rep(-10, 3), + upper = rep(10, 3), parallel = "external") +# to directly call the smcSampler, use +SMCout <- smcSampler(bayesianSetup, initialParticles = 100, iterations = 3) +# plotting results +marginalPlot(SMCout, type = "d", singlePanel = T) +# for options, see +?smcSampler - +# as for all MCMCs, you can also run the SMC via the general runMCMC function +SMCout<-runMCMC(bayesianSetup, sampler = "SMC", + settings = list(nrChains =2, initialParticles = 100, iterations = 3)) +plot(SMCout) } \author{ -Florian Hartig +Florian Hartig, Matthias Speich } diff --git a/BayesianTools/tests/testthat/test-SMC.R b/BayesianTools/tests/testthat/test-SMC.R new file mode 100644 index 0000000..369cc14 --- /dev/null +++ b/BayesianTools/tests/testthat/test-SMC.R @@ -0,0 +1,147 @@ +context("Test SMC and utility functions") # Must be at the top of the file (https://stackoverflow.com/questions/50083521/error-in-xmethod-attempt-to-apply-non-function-in-testthat-test-when) + +library(testthat) +library(BayesianTools) +library(Matching) + +test <- "exact" + +#skip_on_cran() + + + +set.seed(123) + + + +test_that("resampling works", { + + weights = log(c(1,2,3)/6) + size = 10000 + + x1 = replicate(size, BayesianTools:::resample(weights)) + dif = table(x1) - exp(weights) * size * 3 + expect_lt(sum(abs(dif)), 300) + + x2 = replicate(size, resample(weights, method = "residual")) + dif = table(x2) - exp(weights) * size * 3 + expect_lt(sum(abs(dif)), 300) + + x3 = replicate(size, resample(weights, method = "systematic")) + dif = table(x3) - exp(weights) * size * 3 + expect_lt(sum(abs(dif)), 300) + +}) + +test_that("adaptive scheme works", { + weights <- log(rep(1/100, 100)) + ess <- 1 / sum(exp(2 * weights)) + posteriorValues <- log(rnorm(100,10,1)) + posteriorValues <- posteriorValues - BayesianTools:::logSumExp(posteriorValues) + importanceValues <- log(runif(100,0,1)) + + ess.factor <- 0.95 + + # Initial iteration: start with uniform importance distribution + curExp <- 0 + interDist <- importanceValues + + inter.out <- beta.search(ess=ess, target.ess=ess.factor * ess, posteriorValues = posteriorValues, importanceValues = importanceValues, oldInter = interDist, curWeights = weights, curExp = curExp, tol=1) + + curExp <- inter.out$newExp + weights <- inter.out$weights + interDist <- inter.out$interDist + new.ess <- inter.out$ess + + # First iteration: tests whether new ESS is less than original one, and within the tolerance limit specified as argument to beta.search (tol = 1) + + expect_lt(new.ess,ess) + expect_gt(new.ess,(ess*ess.factor - 1)) + + inter.out <- beta.search(ess=new.ess, target.ess=ess.factor * ess, posteriorValues = posteriorValues, importanceValues = importanceValues, oldInter = interDist, curWeights = weights, curExp = curExp, tol=1) + new.new.ess <- inter.out$ess + + # Second iteration: tests whether new ESS is less than previous one, and within the tolerance limit specified as argument to beta.search (tol = 1) + + expect_lt(new.new.ess,new.ess) + expect_gt(new.new.ess,(ess*ess.factor^2 - 1)) +}) + +test_that("mutation works", { + + # Create test BayesianSetup + ll = function(x) sum(dnorm(x, log = T)) + setup = createBayesianSetup(ll, lower = c(-10), upper = c(10)) + + particles <- runif(n=5000, min=-10, max=10) + particles <- matrix(particles, ncol=1) + importanceDensity <- function(x){return(dunif(x, min=-10, max=10, log = TRUE))} + + ### Testing *Metropolis* implementation + + posteriorValues <- setup$posterior$density(particles) + # Initialize proposal Generator + proposalGenerator <- createProposalGenerator(covariance=var(particles)) + mutate.test <- mutate(setup = setup, particles = particles, proposalGenerator = proposalGenerator, posteriorValues = posteriorValues, importanceDensity = importanceDensity, method="Metropolis", steps=50, proposalScale = 0.5) + + # Outputs should be normally distributed after a couple of steps + ks.pval <- ks.test(mutate.test$particles, rnorm(n=length(test.smc$particles)))$p.value + expect_gt(ks.pval,0.05) + + ### Testing *Differential Evolution* implementation + + mutate.test <- mutate(setup = setup, particles = particles, proposalGenerator = proposalGenerator, posteriorValues = posteriorValues, importanceDensity = importanceDensity, method="DE", steps=50, proposalScale = 0.5) + + # Outputs should be normally distributed after a couple of steps + ks.pval <- ks.test(mutate.test$particles, rnorm(n=length(test.smc$particles)))$p.value + expect_gt(ks.pval,0.05) +}) + +test_that("univariate SMC sampler works", { + ll = function(x) sum(dnorm(x, log = T)) + setup = createBayesianSetup(ll, lower = c(-10), upper = c(10)) + + initialParticles <- list() + initialParticles$particles <- runif(n=5000, min=-10, max=10) + initialParticles$density <- function(x){return(dunif(x, min=-10, max=10, log = TRUE))} + + # Test with *Metropolis* mutation + test.smc <- smcSampler(setup, initialParticles = initialParticles, mutate.method ="Metropolis",sampling="systematic",adaptive=TRUE, ess.factor=0.99, resamplingSteps = 1) + ks.pval <- ks.boot(rnorm(n=length(test.smc$particles)), test.smc$particles)$ks.boot.pvalue + expect_gt(ks.pval,0.05) + + # Test with *Differential Evolution* mutation + test.smc <- smcSampler(setup, initialParticles = initialParticles, mutate.method ="DE",sampling="systematic",adaptive=TRUE, ess.factor=0.99, resamplingSteps = 1) + ks.pval <- ks.boot(rnorm(n=length(test.smc$particles)), test.smc$particles)$ks.boot.pvalue + expect_gt(ks.pval,0.05) +}) + +test_that("2D sampler works", { + ll = function(x) sum(dnorm(x, log = T)) + setup = createBayesianSetup(ll, lower = c(-10,-10), upper = c(10,10)) + + initialParticles <- list() + part1 <- runif(n=5000, min=-10, max=10) + part2 <- runif(n=5000, min=-10, max=10) + initialParticles$particles <- data.frame(part1,part2) + initialParticles$density <- function(x){return(rep(dunif(1, min=-10, max=10, log = TRUE), nrow(x)))} + + ### Test with *Metropolis* mutation + test.smc <- smcSampler(setup, initialParticles = initialParticles, mutate.method ="Metropolis",sampling="systematic",adaptive=TRUE, ess.factor=0.99, resamplingSteps = 1) + + #ks.pval1 <- ks.test(test.smc$particles[,1], rnorm(n=length(test.smc$particles)))$p.value + ks.pval1 <- ks.boot(rnorm(n=length(test.smc$particles[,1])), test.smc$particles[,1])$ks.boot.pvalue + expect_gt(ks.pval1,0.05) + + ks.pval2 <- ks.boot(rnorm(n=length(test.smc$particles[,2])), test.smc$particles[,2])$ks.boot.pvalue + expect_gt(ks.pval2,0.05) + + ### Test with *Differential Evolution* mutation + test.smc <- smcSampler(setup, initialParticles = initialParticles, mutate.method ="DE",sampling="systematic",adaptive=TRUE, ess.factor=0.99, resamplingSteps = 1) + + ks.pval1 <- ks.boot(rnorm(n=length(test.smc$particles[,1])), test.smc$particles[,1])$ks.boot.pvalue + expect_gt(ks.pval1,0.05) + + ks.pval2 <- ks.boot(rnorm(n=length(test.smc$particles[,2])), test.smc$particles[,2])$ks.boot.pvalue + expect_gt(ks.pval2,0.05) +}) diff --git a/BayesianTools/tests/testthat/test-settingsDefault.R b/BayesianTools/tests/testthat/test-settingsDefault.R index b88081d..04afbb2 100644 --- a/BayesianTools/tests/testthat/test-settingsDefault.R +++ b/BayesianTools/tests/testthat/test-settingsDefault.R @@ -14,8 +14,6 @@ test_that("Default works in principle",{ } ) - - test_that("Wrong inputs are caught",{ skip_on_cran() diff --git a/Development/SMC/SMC-basicTests.R b/Development/SMC/SMC-basicTests.R new file mode 100644 index 0000000..7060437 --- /dev/null +++ b/Development/SMC/SMC-basicTests.R @@ -0,0 +1,53 @@ + +ll = function(x){ + #for(i in 1:10000) mean(rnorm(10)) + sum(dnorm(x, log = T)) +} +setup = createBayesianSetup(ll, lower = c(-10, -10), upper = c(10, 10)) + +settings = list(iterations = 40000) +out = runMCMC(bayesianSetup = setup, settings = settings) +ks.test(getSample(out, start = 5000)[,1], msm::ptnorm, 0, 1 , -10, 10) +plot(out) + + +initialPart = matrix(rnorm(10000), ncol = 2) + + +out = smcSampler(bayesianSetup = setup, initialParticles = initialPart, iterations = 5, adaptive= F, ess.limit = 20000, mutate = "D") + +par(mfrow = c(1,2)) +hist(initialPart, xlim = c(-4,4)) +hist(out$particles, xlim = c(-4,4)) + +sd(initialPart) +sd(out$particles) +1/sqrt(2) + + + +qqnorm(out$particles) +ks.test(out$particles[,1], msm::ptnorm, 0, 1 , -10, 10) + +out$info$diagnostics +out +summary(out) +plot(out) + + +ptm = proc.time() + +setup = createBayesianSetup(ll, lower = c(-10, -10), upper = c(10, 10)) + +out = smcSampler(bayesianSetup = setup, initialParticles = 2000, iterations = 2, ess.limit = 10000000, diagnostics = function(x)sum(mean(x))) + +proc.time() - ptm + + +ptm = proc.time() + +setup = createBayesianSetup(ll, lower = c(-10, -10), upper = c(10, 10), parallel = T) + +out = smcSampler(bayesianSetup = setup, initialParticles = 2000, iterations = 2, ess.limit = 10000000, diagnostics = function(x)sum(mean(x))) + +proc.time() - ptm \ No newline at end of file diff --git a/Examples/PosteriorVolume.R b/Examples/PosteriorVolume.R new file mode 100644 index 0000000..017f5d5 --- /dev/null +++ b/Examples/PosteriorVolume.R @@ -0,0 +1,24 @@ +library(BayesianTools) + + + +#x1 is technically unidentifiable, only change sensitivity to x2 +#x2 has maximum at zero, regardless of x1 +ll <- function(x){ + dens = (- x[2]^2 / exp(x[1])) +} + +x <- seq(-3,3,length.out=100) +y <- seq(-3,3,length.out=100) +z <- outer(x, y, Vectorize(function(x,y) ll(c(x,y)))) + +image(x, y, log(-z), xlab = "par1", ylab = "par2") + + +# fitting this likelihood with BT and flat priors +bayesianSetup <- createBayesianSetup(likelihood = ll, lower = rep(-10, 2), upper = rep(10, 2)) +out <- runMCMC(bayesianSetup = bayesianSetup) + +marginalPlot(out, histogram = F, singlePanel = T) +correlationPlot(out) +MAP(out) diff --git a/Publications/SpeichEtAl-SMCForModelCalibration/7_plot_results/process_smc_output.r b/Publications/SpeichEtAl-SMCForModelCalibration/7_plot_results/process_smc_output.r new file mode 100644 index 0000000..01bdb13 --- /dev/null +++ b/Publications/SpeichEtAl-SMCForModelCalibration/7_plot_results/process_smc_output.r @@ -0,0 +1,317 @@ +process.out <- function(dir, bin.width=20, bin.max=8000){ + + # This function takes in the output of SMC experiments (including the MCMC reference runs) and returns + # a list with the following objects: + # settings.conv: a data.frame containing the settings (SMC parameters), time and performance of the SMC runs + # (only for those runs which reached convergence) + # bins: vector of equal-interval bins used for the display of the MCMC reference curves + # min.ref: vector containing the lower limit of the MCMC "polygon" (lowest Ds value at each time step / bin) + # max.ref: vector containing the upper limit of the MCMC "polygon" (highest Ds value at each time step / bin) + + # Arguments: + # dir: path to directory where + + library(BayesianTools) + library(coda) + + # Get a list of all settings combinations + + runs <- list.files(paste0(dir,"5/"), pattern = ".RData", full.names = FALSE) + + # Read in the outputs (samples) for all settings combinations (there are 5 runs for each combination). + # Obtain Gelman-Rubin convergence diagnostic between the 5 runs of each combination. + + conv <- list() + + for(i in 1:length(runs)){ + print(runs[i]) + outs <- list() + + for(j in 1:5){ + lpath <- paste0(dir, j, "/") + loaded <- load(file=paste0(lpath,runs[i])) + outs[[j]] <- mcmc(get(loaded)$particles) + } + + conv[[i]] <- gelman.diag(outs) # From coda package + } + + # For each settings combination, get the maximum Gelman-Rubin score for point estimate and upper CI, as well as multivariate Gelman-Rubin. + + max.gr.point <- max.gr.uci <- gr.multi <- vector("numeric", length(runs)) + + for(i in 1:length(runs)){ + max.gr.point[i] <- max(conv[[i]]$psrf[,1]) + max.gr.uci[i] <- max(conv[[i]]$psrf[,2]) + gr.multi[i] <- conv[[i]]$mpsrf + } + + settings <- expand.grid(mutateSteps = c(10, 2, 20, 30, 5), proposalScale = c(0.01, 0.1, 0.333, 0.5), ess.limit = c(0.5, 0.75, 0.9), particles = c(100000, 20000, 50000, 5000)) + + settings$max.gr.point <- max.gr.point + settings$max.gr.uci <- max.gr.uci + settings$gr.multi <- gr.multi + settings$ind <- 1:nrow(settings) + + + # Keep only those settings combinations where the runs have converged + + settings.conv <- settings[settings$max.gr.point < 1.05 & settings$max.gr.uci < 1.05 & settings$gr.multi < 1.2, ] + + # Read in MCMC reference results and SMC results + mcmc.out.np1 <- read.table(paste0(dir,"1/mcmcRef_np.txt"), header=TRUE) + mcmc.out.np2 <- read.table(paste0(dir,"2/mcmcRef_np.txt"), header=TRUE) + mcmc.out.np3 <- read.table(paste0(dir,"3/mcmcRef_np.txt"), header=TRUE) + mcmc.out.np4 <- read.table(paste0(dir,"4/mcmcRef_np.txt"), header=TRUE) + mcmc.out.np5 <- read.table(paste0(dir,"5/mcmcRef_np.txt"), header=TRUE) + + smc.out <- list() + for(j in 1:5){ + smc.out[[j]] <- read.table(paste0(dir, j, "/smc_out.txt"), header=TRUE) + } + + avg.time <- min.time <- max.time <- avg.d <- min.d <- max.d <- vector("numeric", nrow(settings.conv)) + + + + for(i in 1:nrow(settings.conv)){ + time.vec <- d.vec <- vector("numeric", 5) + ind <- which(smc.out[[1]]$particles == settings.conv[i,]$particles & smc.out[[1]]$ess.limit == settings.conv[i,]$ess.limit & smc.out[[1]]$proposalScale == settings.conv[i,]$proposalScale & smc.out[[1]]$mcmcSteps == settings.conv[i,]$mutateSteps) + for(j in 1:5){ + print(c(settings.conv[i,]$particles, settings.conv[i,]$ess.limit, settings.conv[i,]$proposalScale, settings.conv[i,]$mutateSteps)) + time.vec[j] <- smc.out[[j]]$time[ind] + d.vec[j] <- smc.out[[j]]$d[ind] + } + avg.time[i] <- mean(time.vec) + min.time[i] <- min(time.vec) + max.time[i] <- max(time.vec) + avg.d[i] <- mean(d.vec) + min.d[i] <- min(d.vec) + max.d[i] <- max(d.vec) + } + + settings.conv$avg.time <- avg.time + settings.conv$min.time <- min.time + settings.conv$max.time <- max.time + settings.conv$avg.d <- avg.d + settings.conv$min.d <- min.d + settings.conv$max.d <- max.d + + ################### + # Aggregate the 5 MCMC reference runs + bins <- seq(0,bin.max,by=bin.width) + min.ref <- max.ref <- vector("numeric", (length(bins)-1)) + + for(i in 1:length(min.ref)){ + + all.inbin <- c(mcmc.out.np1$distance[mcmc.out.np1$time > bins[i] & mcmc.out.np1$time <= bins[i+1]], + mcmc.out.np2$distance[mcmc.out.np2$time > bins[i] & mcmc.out.np2$time <= bins[i+1]], + mcmc.out.np3$distance[mcmc.out.np3$time > bins[i] & mcmc.out.np3$time <= bins[i+1]], + mcmc.out.np4$distance[mcmc.out.np4$time > bins[i] & mcmc.out.np4$time <= bins[i+1]], + mcmc.out.np5$distance[mcmc.out.np5$time > bins[i] & mcmc.out.np5$time <= bins[i+1]]) + + min.ref[i] <- min(all.inbin) + max.ref[i] <- max(all.inbin) + } + + out <- list(settings.conv, bins, min.ref, max.ref) + return(out) +} + +################################################################# +################################################################# + +plot.smc.out <- function(out.list, xlim = c(0,7000), ylim = c(0,0.2), main="Model", bars=TRUE, symbol.scale = 1, xlog = "no", polygon.col = adjustcolor("blue", alpha.f = 0.2)){ + +# This function plots one panel of Fig. 3 in the manuscript of Speich et al. + + settings.conv <- out.list[[1]] + bins <- out.list[[2]] + min.ref <- out.list[[3]] + max.ref <- out.list[[4]] + + + settings.conv$col <- NA + settings.conv$col[settings.conv$proposalScale==0.01] <- "yellow" + settings.conv$col[settings.conv$proposalScale==0.1] <- "orange" + settings.conv$col[settings.conv$proposalScale==0.333] <- "red" + settings.conv$col[settings.conv$proposalScale==0.5] <- "black" + + + settings.conv$cex <- NA + settings.conv$cex[settings.conv$particles==5000] <- 0.25 + settings.conv$cex[settings.conv$particles==20000] <- 0.5 + settings.conv$cex[settings.conv$particles==50000] <- 1 + settings.conv$cex[settings.conv$particles==100000] <- 2 + settings.conv$cex <- settings.conv$cex * symbol.scale + + settings.conv$fill <- NA + settings.conv$fill[settings.conv$ess.limit==0.5] <- "#edf8fb" + settings.conv$fill[settings.conv$ess.limit==0.75] <- "#b2e2e2" + settings.conv$fill[settings.conv$ess.limit==0.9] <- "#66c2a4" + #settings.conv$fill[settings.conv$ess.factor==0.99] <- "#238b45" + + settings.conv$pch <- NA + settings.conv$pch[settings.conv$mutateSteps == 2] <- 21 + settings.conv$pch[settings.conv$mutateSteps == 5] <- 22 + settings.conv$pch[settings.conv$mutateSteps == 10] <- 23 + settings.conv$pch[settings.conv$mutateSteps == 20] <- 24 + settings.conv$pch[settings.conv$mutateSteps == 30] <- 25 + + if(xlog=="ln"){ + avg.time <- log(settings.conv$avg.time) + if(bins[1]==0) bins[1] <- 1 + bins <- log(bins) + min.time <- log(settings.conv$min.time) + max.time <- log(settings.conv$max.time) + if(xlim[2] >= 10000){ + x.labels <- seq(round(xlim[1],-3), xlim[2], by=2000) + x.ticks <- log(seq(round(xlim[1],-2), xlim[2], by=200)) + } else { + x.labels <- seq(round(xlim[1],-3), xlim[2], by=1000) + x.ticks <- log(seq(round(xlim[1],-2), xlim[2], by=100)) + } + x.pos <- log(x.labels) + xlim <- log(xlim) + } else if(xlog=="log10"){ + avg.time <- log10(settings.conv$avg.time) + if(bins[1]==0) bins[1] <- 1 + bins <- log10(bins) + min.time <- log10(settings.conv$min.time) + max.time <- log10(settings.conv$max.time) + if(xlim[2] >= 10000){ + x.labels <- seq(round(xlim[1],-3), xlim[2], by=2000) + x.ticks <- log10(seq(round(xlim[1],-2), xlim[2], by=200)) + } else { + x.labels <- seq(round(xlim[1],-3), xlim[2], by=1000) + x.ticks <- log10(seq(round(xlim[1],-2), xlim[2], by=100)) + } + x.pos <- log10(x.labels) + xlim <- log10(xlim) + } else if(xlog=="no"){ + avg.time <- settings.conv$avg.time + min.time <- settings.conv$min.time + max.time <- settings.conv$max.time + if(xlim[2] >= 20000){ + x.labels <- seq(round(xlim[1],-3), xlim[2], by=2000) + x.ticks <- seq(round(xlim[1],-2), xlim[2], by=200) + } else { + x.labels <- seq(round(xlim[1],-3), xlim[2], by=1000) + x.ticks <- seq(round(xlim[1],-2), xlim[2], by=200) + } + x.pos <- x.labels + } + + plot(avg.time, settings.conv$avg.d, type="n", ylim=ylim, xlim=xlim, xlab="Runtime [s]", ylab="Distance from reference", main = main, las=1, axes=FALSE, cex.main=0.8) + #polygon.col <- adjustcolor("blue", alpha.f = 0.2) + polygon(x = c(bins[1:(length(bins)-1)],rev(bins[1:(length(bins)-1)])), y=c(min.ref, rev(max.ref)), col=polygon.col, border=NA) + if(bars){ + segments(x0=min.time, y0=settings.conv$avg.d, x1=max.time, lwd=0.5) + segments(x0 = avg.time, y0 = settings.conv$min.d, y1 = settings.conv$max.d, lwd=0.5) + } + points(avg.time, settings.conv$avg.d, pch=settings.conv$pch, col=settings.conv$col, cex=settings.conv$cex, bg=settings.conv$fill) + + axis(side = 1, at=x.pos, labels = x.labels, las=3, cex.axis = 0.8) + axis(side=1, lwd = 0, at=x.ticks, labels = rep("", length(x.ticks)), tck = -0.02, lwd.ticks = 0.5) + axis(side = 2, las = 1, cex.axis = 0.8) +} + + +processed.VSEMa <- process.out("/Users/mspeich/Documents/bayes/output/VSEMa") +processed.VSEMb <- process.out("/Users/mspeich/Documents/bayes/output/VSEMb") +processed.threePGN <- process.out("/Users/mspeich/Documents/bayes/output/threePGN") +processed.threePGN_sleep <- process.out("/Users/mspeich/Documents/bayes/output/threePGN_sleep", bin.max=30000, bin.width=300) + +proc.out <- list(processed.VSEMa, processed.VSEMb, processed.threePGN, processed.threePGN_sleep) +save(proc.out, file="/Users/mspeich/Documents/bayes/output/processed.RData") + +plot.smc.out(processed.VSEMa, main="VSEM - Strong correlation", xlim=c(10,2000), bars=TRUE, symbol.scale = 1.5, xlog="ln") +plot.smc.out(processed.VSEMa, main="VSEM - Strong correlation", xlim=c(10,2000), bars=TRUE, symbol.scale = 1.5, xlog="log10") +plot.smc.out(processed.VSEMb, main="VSEM - No strong correlation", xlim=c(10,5000), bars=FALSE, symbol.scale=1.5, xlog="log10") +plot.smc.out(processed.threePGN, main="3-PGN", xlim=c(200,10000), bars=TRUE, symbol.scale = 1.5, xlog = "log10") + +processed.threePGN_sleep <- process.out("/Users/mspeich/Documents/bayes/output/threePGN_sleep") + + +################################### +### The following code generates the figure shown as Fig. 3 in Speich et al. + +pdf(file ="/Users/mspeich/Documents/bayes/figs/smc_compare_results_linTime.pdf", height = 8, width = 9) +#png(file ="/Users/mspeich/Documents/bayes/figs/smc_compare_results_linTime.png", height = 8, width = 9,units="in") + +#par(mfrow=c(2,2)) +par(mar=c(3.1, 3.1, 1.1, 1.1)) +par(fig=c(0.05, 0.525, 0.6, 1)) +#plot.smc.out(processed.VSEMa, main="VSEM - Strong correlation", xlim=c(10,5000), bars=FALSE, symbol.scale = 1, xlog="log10", polygon.col = "lightblue") +plot.smc.out(processed.VSEMa, main="VSEM - Strong correlation", xlim=c(10,5000), bars=FALSE, symbol.scale = 1, xlog="no", polygon.col = "lightblue") + +par(fig=c(0.525, 1, 0.6, 1), new=TRUE) +#plot.smc.out(processed.VSEMb, main="VSEM - No strong correlation", xlim=c(10,5000), bars=FALSE, symbol.scale=1, xlog="log10", polygon.col = "lightblue") +plot.smc.out(processed.VSEMb, main="VSEM - No strong correlation", xlim=c(10,5000), bars=FALSE, symbol.scale=1, xlog="no", polygon.col = "lightblue") + +par(fig=c(0.05, 0.525, 0.2, 0.6), new=TRUE) +#plot.smc.out(processed.threePGN, main="3-PGN", xlim=c(20,10000), bars=FALSE, symbol.scale = 1, xlog = "log10", polygon.col = "lightblue") +plot.smc.out(processed.threePGN, main="3-PGN", xlim=c(20,8000), bars=FALSE, symbol.scale = 1, xlog = "no", polygon.col = "lightblue") + + +par(fig=c(0.525, 1, 0.2, 0.6), new=TRUE) +#plot.smc.out(processed.sleep, main="3-PGN sleep", xlim=c(10,20000), bars=FALSE, symbol.scale = 1, xlog = "log10", polygon.col = "lightblue") # TEMP +plot.smc.out(processed.threePGN_sleep, main="3-PGN sleep", xlim=c(10,24000), bars=FALSE, symbol.scale = 1, xlog = "no", polygon.col = "lightblue") # TEMP + +## LEGEND +par(fig=c(0, 1, 0, 0.15),new=TRUE) +par(mar=c(0,0,0,0)) +plot(1,1, type = "n", axes=FALSE, xlab="", ylab="", xlim = c(0,7000), ylim=c(-0.01, 0.05)) + +points(x=0, y=0.03, pch=21, cex=0.25, col="black", bg="white") +points(x=0, y=0.02, pch=21, cex=0.5, col="black", bg="white") +points(x=0, y=0.01, pch=21, cex=1, col="black", bg="white") +points(x=0, y=0, pch=21, cex=2, col="black", bg="white") +text(x=200, y=0.04, labels="N", font=2) +text(x=400, y=0.03, labels="5000", adj=0, cex=0.8) +text(x=400, y=0.02, labels="20000", adj=0, cex=0.8) +text(x=400, y=0.01, labels="50000", adj=0, cex=0.8) +text(x=400, y=0, labels="100000", adj=0, cex=0.8) + +points(x=1500, y=0.03, pch=21, col="black", bg="#edf8fb") +points(x=1500, y=0.015, pch=21, col="black", bg="#b2e2e2") +points(x=1500, y=0, pch=21, col="black", bg="#66c2a4") +text(x=1700, y=0.04, labels="a", font=2) +text(x=1900, y=0.03, labels = "0.5", adj=0, cex = 0.8) +text(x=1900, y=0.015, labels = "0.75", adj=0, cex = 0.8) +text(x=1900, y=0, labels = "0.9", adj=0, cex = 0.8) + +points(x=3000, y = 0.03, pch = 21, col = "black", bg = "white") +points(x=3000, y = 0.03 * 0.75, pch = 22, col = "black", bg = "white") +points(x=3000, y = 0.03 * 0.5, pch = 23, col = "black", bg = "white") +points(x=3000, y = 0.03 * 0.25, pch = 24, col = "black", bg = "white") +points(x=3000, y = 0, pch = 25, col = "black", bg = "white") +text(x=3200, y=0.04, labels = "S", font = 2) +text(x=3400, y = 0.03, labels = "2", adj = 0, cex = 0.8) +text(x=3400, y = 0.03 * 0.75, labels = "5", adj = 0, cex = 0.8) +text(x=3400, y = 0.03 * 0.5, labels = "10", adj = 0, cex = 0.8) +text(x=3400, y = 0.03 * 0.25, labels = "20", adj = 0, cex = 0.8) +text(x=3400, y = 0, labels = "30", adj = 0, cex = 0.8) + +points(x = 4500, y = 0.03, pch = 21, col = "yellow", bg = "white") +points(x = 4500, y = 0.02, pch = 21, col = "orange", bg = "white") +points(x = 4500, y = 0.01, pch = 21, col = "red", bg = "white") +points(x = 4500, y = 0.00, pch = 21, col = "black", bg = "white") +text(x = 4700, y = 0.04, labels = expression(gamma), font = 2) +text(x = 4900, y = 0.03, labels = "0.01", adj = 0, cex = 0.8) +text(x = 4900, y = 0.02, labels = "0.1", adj = 0, cex = 0.8) +text(x = 4900, y = 0.01, labels = "0.333", adj = 0, cex = 0.8) +text(x = 4900, y = 0, labels = "0.5", adj = 0, cex = 0.8) + +rect(xleft=6000, xright = 7000, ybottom = 0.02, ytop=0.03, border = NA, col = "lightblue") +text(x = 6000, y = 0.015, labels = "MCMC reference", cex = 0.8, adj = 0) + + +## Annotations +par(fig=c(0,1,0,1), new=TRUE) +par(mar=c(0,0,0,0)) +plot(1, 1, type="n", axes = FALSE, xlab="", ylab="", xlim=c(0,1), ylim=c(0,1)) +text(x = 0.52, y = 0.15, labels = "Runtime [s]") +text(x = 0.01, y = 0.6, labels = "Distance from reference [-]", srt = 90) + +dev.off() \ No newline at end of file