From 6ccf229a4ab27dacc4c816ac865229825267b8d7 Mon Sep 17 00:00:00 2001 From: Florian Hartig Date: Fri, 13 Jul 2018 16:08:04 +0200 Subject: [PATCH 01/42] distance function --- BayesianTools/NAMESPACE | 1 + {Test => BayesianTools/R}/Distance.R | 19 ++------- .../inst/examples/correlationPlotHelp.R | 2 +- .../inst/examples/getSampleDistanceHelp.R | 14 +++++++ .../inst/examples/marginalPlotHelp.R | 2 +- BayesianTools/inst/examples/tracePlotHelp.R | 2 +- BayesianTools/man/getSampleDistance.Rd | 41 +++++++++++++++++++ Examples/PosteriorVolume.R | 24 +++++++++++ 8 files changed, 86 insertions(+), 19 deletions(-) rename {Test => BayesianTools/R}/Distance.R (86%) create mode 100644 BayesianTools/inst/examples/getSampleDistanceHelp.R create mode 100644 BayesianTools/man/getSampleDistance.Rd create mode 100644 Examples/PosteriorVolume.R diff --git a/BayesianTools/NAMESPACE b/BayesianTools/NAMESPACE index 0681a97..6776f24 100644 --- a/BayesianTools/NAMESPACE +++ b/BayesianTools/NAMESPACE @@ -66,6 +66,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 86% rename from Test/Distance.R rename to BayesianTools/R/Distance.R index 4f06e8d..d513488 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"){ @@ -72,23 +75,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/inst/examples/correlationPlotHelp.R b/BayesianTools/inst/examples/correlationPlotHelp.R index 9522ea6..aa64e23 100644 --- a/BayesianTools/inst/examples/correlationPlotHelp.R +++ b/BayesianTools/inst/examples/correlationPlotHelp.R @@ -22,4 +22,4 @@ correlationPlot(out, scaleCorText = FALSE, method = "kendall") correlationPlot(out, scaleCorText = FALSE, method = "kendall", whichParameters = c(1,2)) ## additional parameters can be passed to getSample (see ?getSample for further information) -correlationPlot(out, scaleCorText = FALSE, thin = 100, start = 200) \ No newline at end of file +correlationPlot(out, scaleCorText = FALSE, thin = 100, start = 200) diff --git a/BayesianTools/inst/examples/getSampleDistanceHelp.R b/BayesianTools/inst/examples/getSampleDistanceHelp.R new file mode 100644 index 0000000..b78c721 --- /dev/null +++ b/BayesianTools/inst/examples/getSampleDistanceHelp.R @@ -0,0 +1,14 @@ +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") diff --git a/BayesianTools/inst/examples/marginalPlotHelp.R b/BayesianTools/inst/examples/marginalPlotHelp.R index 6366a7d..fdb6b90 100644 --- a/BayesianTools/inst/examples/marginalPlotHelp.R +++ b/BayesianTools/inst/examples/marginalPlotHelp.R @@ -28,4 +28,4 @@ marginalPlot(out, plotPrior = FALSE) marginalPlot(out, singlePanel = TRUE, scale=TRUE, col = c("red", "blue"), thin = 100) marginalPlot(out, singlePanel = TRUE, scale=TRUE, - col = c("red", "blue"), thin = 100, start = 500) \ No newline at end of file + col = c("red", "blue"), thin = 100, start = 500) diff --git a/BayesianTools/inst/examples/tracePlotHelp.R b/BayesianTools/inst/examples/tracePlotHelp.R index 127e441..cc50aa4 100644 --- a/BayesianTools/inst/examples/tracePlotHelp.R +++ b/BayesianTools/inst/examples/tracePlotHelp.R @@ -11,4 +11,4 @@ tracePlot(sampler = out, thin = 50) # additional parameters can be passed on to getSample (see help) tracePlot(sampler = out, thin = 10, start = 500) # select parameter by index -tracePlot(sampler = out, thin = 10, start = 500, whichParameters = 2) \ No newline at end of file +tracePlot(sampler = out, thin = 10, start = 500, whichParameters = 2) diff --git a/BayesianTools/man/getSampleDistance.Rd b/BayesianTools/man/getSampleDistance.Rd new file mode 100644 index 0000000..3d3d1fd --- /dev/null +++ b/BayesianTools/man/getSampleDistance.Rd @@ -0,0 +1,41 @@ +% 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") +} 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) From b17e112c33604582a56c262e9a0840745f12db43 Mon Sep 17 00:00:00 2001 From: Florian Hartig Date: Tue, 31 Jul 2018 16:19:33 +0200 Subject: [PATCH 02/42] updates to SMC --- BayesianTools/R/SMC.R | 31 +++++++++++++++++++++++-------- 1 file changed, 23 insertions(+), 8 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index 10fd53b..b1d8c09 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -49,21 +49,33 @@ smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, for (i in 1:iterations){ + # 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 + posterior = setup$posterior$density(particles, returnAll = T) likelihoodValues <- posterior[,2] + + ll = likelihoodValues - max(likelihoodValues, na.rm = T) - # idea - adjust (1/iterations) such that always approx 30% of particles are maintain - #level = sort(likelihoodValues)[acceptanceTarget] - #best = likelihoodValues + llCutoff = sort(ll)[acceptanceTarget] - ll = likelihoodValues - max(likelihoodValues, na.rm = T) + relativeL = exp(likelihoodValues - max(likelihoodValues, na.rm = T))^(1/iterations) - llCutoff = sort(ll)[acceptanceTarget] + # Re?sampling step + # Option 1) Always resample (default BT) + # Option 2) Keep weights and resample if ESS < something - relativeL = exp(likelihoodValues - max(likelihoodValues, na.rm = T))^(1/iterations) + # Question: last step complete resampling, so that all particles are the same? sel = sample.int(n=length(likelihoodValues), size = length(likelihoodValues), replace = T, prob = relativeL) @@ -73,6 +85,11 @@ smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, if (numPar == 1) particles = matrix(particles, ncol = 1) + # Mutate + + # 1) When to do the resampling + # 2) HOW - currently adaptive Metropolis, could also do DEzs step + if (resampling == T){ if (adaptive == T){ @@ -89,8 +106,6 @@ smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, rejectionRate = rejectionRate + sum(accepted) particles[accepted, ] = particlesProposals[accepted, ] - - } } } From 7c021d4c8663b3ac1dd913d737a6ae5567f5205c Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Wed, 1 Aug 2018 11:11:34 +0200 Subject: [PATCH 03/42] Specified numerical tolerance parameter to avoid crashes It has occurred that getSampleDistance crashed with option "BH" because of very small values in the inverse covariance matrix (see https://stackoverflow.com/questions/22134398/mahalonobis-distance-in-r-error-system-is-computationally-singular). The argument tol determines which (very small) values should be considered equal to zero. --- BayesianTools/R/Distance.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/BayesianTools/R/Distance.R b/BayesianTools/R/Distance.R index d513488..c644e82 100644 --- a/BayesianTools/R/Distance.R +++ b/BayesianTools/R/Distance.R @@ -30,7 +30,7 @@ getSampleDistance <- function(sample1, sample2, type = "KL"){ # The following code is copied from package fpc aggregatesigma <- (Sigma1+Sigma2)/2 - d1 <- mahalanobis(mu1,mu2,aggregatesigma)/8 + 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 d2 <- log(det(as.matrix(aggregatesigma))/sqrt(det(as.matrix(Sigma1))* det(as.matrix(Sigma2))))/2 out <- d1+d2 From 0c34c9b127c97a45f60eeb7bd4b1244a7626ddaa Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Wed, 1 Aug 2018 11:56:24 +0200 Subject: [PATCH 04/42] Modify particle resampling scheme Several modifications to the particle resampling step: - generate a sequence of intermediary distributions between prior and posterior - use weights attributed to particles - calculate Effective Sample Size (ESS) to determine when particles should be resampled - dynamically adapt the number of iterations --- BayesianTools/R/SMC.R | 231 +++++++++++++++++++++++++++++++++++++----- 1 file changed, 205 insertions(+), 26 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index b1d8c09..eb42e11 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -2,28 +2,45 @@ #' @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 reference object of class MCMCsampler, used as a ground truth to assess performance of the SMC sampler #' @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 if new particles should be created at each iteration #' @param resamplingSteps how many resampling (MCMC) steps between the iterations #' @param proposal optional proposal class +#' @param exponents series of exponents to build the intermediate distributions #' @param adaptive 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 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 of effective sample size below which resampling is done. By default, the value is set to half the number of particles. To resample at each step, use a value >= the number of particles. +#' @param lastResample Iteration (starting from the end) at which particle resampling is forced. To deactivate this, set to a value < 0 +#' @param pars.lower Optional: vector contaiing the mimimum values of calibration parameters. Used to initialize the proposal function +#' @param pars.upper Optional: vector contaiing the maximum values of calibration parameters. Used to initialize the proposal function #' @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. #' @export #' @example /inst/examples/SMCHelp.R -smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, resampling = T, resamplingSteps = 2, proposal = NULL, adaptive = T, proposalScale = 0.5){ +smcSampler <- function(bayesianSetup, reference=NULL, initialParticles = 1000, iterations = 10, resampling = T, resamplingSteps = 2, proposal = NULL, exponents = NULL, adaptive = T, proposalScale = 0.5, x=3.11, m=7E-08, sampling="multinomial",ess.limit=NULL,lastResample = 1,pars.lower=NULL, pars.upper=NULL){ if(resamplingSteps < 1) stop("SMC error, resamplingSteps can't be < 1") - - setup <- checkBayesianSetup(bayesianSetup) info = list() - info$resamplingAcceptance = matrix(nrow = iterations, ncol = resamplingSteps) + info$resamplingAcceptance = as.data.frame(matrix(nrow = iterations, ncol = resamplingSteps)) #Using DF instead of matrix because number of iterations may change due to adaptive algorithm info$survivingParticles = rep(NA, iterations) + info$ess.vec <- vector("numeric", iterations) + info$dist.vec <- vector("numeric", iterations) + info$res.of <- rep(FALSE, iterations) + info$res.ess <- rep(FALSE, iterations) + weights <- vector("numeric", initialParticles) + oldweights <- vector("numeric", initialParticles) + oldll <- vector("numeric", initialParticles) + oldInter <- vector("numeric", initialParticles) + expTrack <- vector("numeric", initialParticles) + + setup <- checkBayesianSetup(bayesianSetup) if(class(initialParticles) == "numeric"){ initialParticles = bayesianSetup$prior$sampler(initialParticles) @@ -32,25 +49,51 @@ smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, if (any(is.infinite(setup$prior$density(initialParticles)))) stop("initialParticles outside prior range") particles <- initialParticles - rejectionRate = 0 - particleSize = nrow(initialParticles) - - acceptanceTarget = round(particleSize / 2) + if(is.null(ess.limit)){ess.limit <- round(particleSize) * 0.5} posterior = matrix(nrow = particleSize, ncol = 3) - numPar <- ncol(initialParticles) if (is.null(proposal)) proposalGenerator = createProposalGenerator(rep(40,numPar)) usedUp = 0 - for (i in 1:iterations){ + # 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] + + icount <- 1 # Iterations counter + + # All particles are given equal weight at initialization + weights[1:length(weights)] <- log(1/particleSize) + + 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, message = TRUE)} + + #for (i in 1:iterations){ + while(icount <= length(exponents)){ + # Using a while loop instead of for because in the adaptive algorithm, the number of iterations may increase # Reweighting - + # idea - adjust (1/iterations) such that always approx 30% of particles are maintain # level = sort(likelihoodValues)[acceptanceTarget] # best = likelihoodValues @@ -60,28 +103,116 @@ smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, # https://github.com/florianhartig/BayesianTools/issues/23 posterior = setup$posterior$density(particles, returnAll = T) + priorValues <- posterior[,3] + likelihoodValues <- posterior[, 2] + postValues <- posterior[,1] - likelihoodValues <- posterior[,2] - - ll = likelihoodValues - max(likelihoodValues, na.rm = T) + print(c("weights", head(weights))) + + ## Update weights + # Create intermediary distribution + curExp <- exponents[icount] + interDist <- curExp * postValues + (1-curExp) * priorValues - llCutoff = sort(ll)[acceptanceTarget] + # Calculate new weights + oldweights <- weights + weights <- weights + (interDist - oldInter) - relativeL = exp(likelihoodValues - max(likelihoodValues, na.rm = T))^(1/iterations) + # Normalize (log-)weights so that the sum (of non-logs) equals 1 + sumw <- sum(sort(exp(weights))) + weights <- weights - log(sumw) + + # It may occur that the difference between the new and old weights of a particle is so great that an overflow occurs. + # In this case, the algorithm reduces the difference between the current and next intermediary distribution (i.e. the difference + # between the current and next value of the exponent). As a result, one more iteration is added. + # This is the same that happens when the effective sample size becomes *very* low. + #if(is.infinite(sum(sort(exp(weights))))){ + if(any(is.infinite(weights))){ + print("Infinite weights sum; retrying with new intermediary distribution") + expdiff <- (exponents[icount] - exponents[(icount-1)]) * 0.5 + newExp <- exponents[(icount-1)] + expdiff + exponents <- append(exponents,newExp,after=(icount-1)) + weights <- oldweights + info$res.of[icount] <- TRUE + next + } + + # posterior = setup$posterior$density(particles, returnAll = T) + # likelihoodValues <- posterior[,2] + # ll = likelihoodValues - max(likelihoodValues, na.rm = T) + # llCutoff = sort(ll)[acceptanceTarget] + # relativeL = exp(likelihoodValues - max(likelihoodValues, na.rm = T))^(1/iterations) # Re?sampling step # Option 1) Always resample (default BT) # Option 2) Keep weights and resample if ESS < something + # => for now, option 1 is activated if the argument ess.limit is set to a value >= the number of particles # Question: last step complete resampling, so that all particles are the same? - - sel = sample.int(n=length(likelihoodValues), size = length(likelihoodValues), replace = T, prob = relativeL) - info$survivingParticles[i] = length(unique(sel)) + #sel = sample.int(n=length(likelihoodValues), size = length(likelihoodValues), replace = T, prob = relativeL) + + # Calculate effective sample size (ESS) + ess <- 1 / sum(exp(2 * weights)) + + # Save intermediate distribution for the next iteration + oldInter <- interDist + + if(is.infinite(ess)){print(weights)} - particles = particles[sel,] + # 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 | icount == (length(exponents) - lastResample)){ + oldExp <- curExp + if(sampling == "multinomial"){ + sel <- sample.int(n=length(likelihoodValues), size = length(likelihoodValues), replace = T, prob = exp(weights)) + } else if(sampling == "residual"){ + sel <- residualResampling(weights) + } else if(sampling == "systematic"){ + sel <- systematicResampling(weights) + } else{ + stop("Invalid string for resampling argument") + } + + particles = particles[sel,] + info$survivingParticles[i] = 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 + sumw <- sum(sort(exp(weights))) + weights <- weights - log(sumw) + + # Get posterior values and distribution for new particles + posterior = setup$posterior$density(particles, returnAll = T) + priorValues <- posterior[,3] + likelihoodValues <- posterior[, 2] + postValues <- posterior[,1] + oldInter <- oldExp * postValues + (1-oldExp) * priorValues + + # Determine if exponents need to be modified + if(ess < estar){ + # If yes: + # - discard newly sampled population of particles + # - do not increase counter + # - modify x + + expdiff <- (exponents[icount] - exponents[(icount-1)]) * 0.5 + newExp <- exponents[(icount-1)] + expdiff + exponents <- append(exponents,newExp,after=(icount-1)) + info$res.ess[icount] <- TRUE + next + } + } + + info$ess.vec[icount] <- ess + + icount <- icount + 1 + curExp <- m * (icount*200/iterations)^x + oldExp <- curExp if (numPar == 1) particles = matrix(particles, ncol = 1) @@ -99,19 +230,20 @@ smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, 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)) + #jumpProb <- exp(setup$posterior$density(particlesProposals) - likelihoodValues[sel])^(i/iterations) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) + jumpProb <- exp(setup$posterior$density(particlesProposals) - likelihoodValues) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) accepted <- jumpProb > runif(length(jumpProb), 0 ,1) - rejectionRate = rejectionRate + sum(accepted) - - particles[accepted, ] = particlesProposals[accepted, ] + particles[accepted, ] = particlesProposals[accepted, ] + info$resamplingAcceptance[(icount-1),j] <- sum(accepted)/particleSize } } + info$dist.vec[(icount-1)] <- getSampleDistance(particles, getSample(reference),type="BH") } info$rejectionRate = rejectionRate / (iterations * resamplingSteps) - + out = list( setup = setup, initialParticles = initialParticles, @@ -125,3 +257,50 @@ smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, return(out) } + + +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){ + #print(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)]) + print(cumu.weights) + 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 + print(c(i, U[i])) + + new.parts.sort[i] <- tail(which(cumu.weights <= U[i]),1) + print(c("new.parts.sort", new.parts.sort[i])) + # Match sorted ranks back to original order of weights + new.parts[i] <- which(rank.weights==new.parts.sort[i]) + } + + return(new.parts) +} \ No newline at end of file From 83ef7b52851a733d49632791a70864201f022ef0 Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Tue, 7 Aug 2018 19:20:31 +0200 Subject: [PATCH 05/42] Add basic Differential Evolution option for SMC --- BayesianTools/R/SMC.R | 50 ++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 47 insertions(+), 3 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index eb42e11..25145a9 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -22,7 +22,7 @@ #' @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, reference=NULL, initialParticles = 1000, iterations = 10, resampling = T, resamplingSteps = 2, proposal = NULL, exponents = NULL, adaptive = T, proposalScale = 0.5, x=3.11, m=7E-08, sampling="multinomial",ess.limit=NULL,lastResample = 1,pars.lower=NULL, pars.upper=NULL){ +smcSampler <- function(bayesianSetup, reference=NULL, initialParticles = 1000, iterations = 10, resampling = T, resamplingSteps = 2, proposal = NULL, exponents = NULL, adaptive = T, proposalScale = 0.5, x=3.11, m=7E-08, sampling="multinomial",ess.limit=NULL,lastResample = 1,pars.lower=NULL, pars.upper=NULL,mutate="Metropolis",b=1e-04){ if(resamplingSteps < 1) stop("SMC error, resamplingSteps can't be < 1") @@ -56,7 +56,9 @@ smcSampler <- function(bayesianSetup, reference=NULL, initialParticles = 1000, i posterior = matrix(nrow = particleSize, ncol = 3) numPar <- ncol(initialParticles) - if (is.null(proposal)) proposalGenerator = createProposalGenerator(rep(40,numPar)) + if(mutate %in% c("Metropolis", "adaptive")){ + if (is.null(proposal)) proposalGenerator = createProposalGenerator(rep(40,numPar)) + } usedUp = 0 @@ -221,7 +223,12 @@ smcSampler <- function(bayesianSetup, reference=NULL, initialParticles = 1000, i # 1) When to do the resampling # 2) HOW - currently adaptive Metropolis, could also do DEzs step - if (resampling == T){ + # 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 + + #if (resampling == T){ + if(mutate=="Metropolis"){ if (adaptive == T){ proposalGenerator = updateProposalGenerator(proposalGenerator, particles) @@ -238,6 +245,43 @@ smcSampler <- function(bayesianSetup, reference=NULL, initialParticles = 1000, i particles[accepted, ] = particlesProposals[accepted, ] info$resamplingAcceptance[(icount-1),j] <- sum(accepted)/particleSize } + } else if(mutate=="DE"){ + # Differential evolution + # For now (test): basic DE algorithm + + for(j in 1:resamplingSteps){ + + # Loop over particles + for(part in 1:particleSize){ + particlesOld <- particles + particleDiff <- rep(0,numPar) + newParts <- rep(0L,2) + + # Sample 2 other particles + while(all(particleDiff==0) | part %in% newParts){ + # 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=particleSize,size=2) + particleDiff <- particlesOld[newParts[1],] - particlesOld[newParts[2],] + } + + particleProp <- particlesOld[part,] + (particleDiff * proposalScale) + runif(numPar,-b,b) + jumpProb <- exp(setup$posterior$density(particleProp) - likelihoodValues[part]) * exp(setup$prior$density(particleProp) - setup$prior$density(particles[part,])) + + if(jumpProb > runif(length(jumpProb), 0 ,1)){ + particles[part,] <- particleProp + if(is.na(info$resamplingAcceptance[(icount-1),j])){ + info$resamplingAcceptance[(icount-1),j] <- 1/particleSize + } else{ + info$resamplingAcceptance[(icount-1),j] <- info$resamplingAcceptance[(icount-1),j] + 1/particleSize + } + } else{ + particles[part,] <- particlesOld[part,] + } + } + } } info$dist.vec[(icount-1)] <- getSampleDistance(particles, getSample(reference),type="BH") } From 85130669d668472554979165b109d582c92600a2 Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Tue, 7 Aug 2018 19:38:12 +0200 Subject: [PATCH 06/42] Fix updating of surviving particles --- BayesianTools/R/SMC.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index 25145a9..d160819 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -180,7 +180,7 @@ smcSampler <- function(bayesianSetup, reference=NULL, initialParticles = 1000, i } particles = particles[sel,] - info$survivingParticles[i] = length(unique(sel)) + info$survivingParticles[icount] = length(unique(sel)) # Set all weights equal weights[1:length(weights)] <- log(1/particleSize) From 6fe7239616caed379534a39ff6d5718beb4870a0 Mon Sep 17 00:00:00 2001 From: Florian Hartig Date: Thu, 9 Aug 2018 12:41:43 +0200 Subject: [PATCH 07/42] SMC tests --- Development/SMC/SMC-basicTests.R | 14 ++++++++++++++ 1 file changed, 14 insertions(+) create mode 100644 Development/SMC/SMC-basicTests.R diff --git a/Development/SMC/SMC-basicTests.R b/Development/SMC/SMC-basicTests.R new file mode 100644 index 0000000..24d808f --- /dev/null +++ b/Development/SMC/SMC-basicTests.R @@ -0,0 +1,14 @@ + +ll = function(x) 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) + + +out = smcSampler(bayesianSetup = setup, initialParticles = 2000, iterations = 1, ess.limit = 10000000) +hist(out$particles) +qqnorm(out$particles) +ks.test(out$particles[,1], msm::ptnorm, 0, 1 , -10, 10) + From 668d2d235e9c47ca32f1f484bfd4b2c6785c0958 Mon Sep 17 00:00:00 2001 From: Florian Hartig Date: Thu, 9 Aug 2018 12:42:33 +0200 Subject: [PATCH 08/42] SMC function changes simplification introduced general diagnostics option --- BayesianTools/R/SMC.R | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index d160819..9171740 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -2,7 +2,6 @@ #' @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 reference object of class MCMCsampler, used as a ground truth to assess performance of the SMC sampler #' @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 if new particles should be created at each iteration @@ -18,11 +17,13 @@ #' @param lastResample Iteration (starting from the end) at which particle resampling is forced. To deactivate this, set to a value < 0 #' @param pars.lower Optional: vector contaiing the mimimum values of calibration parameters. Used to initialize the proposal function #' @param pars.upper Optional: vector contaiing the maximum values of calibration parameters. Used to initialize the proposal function +#' @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. #' @export #' @example /inst/examples/SMCHelp.R -smcSampler <- function(bayesianSetup, reference=NULL, initialParticles = 1000, iterations = 10, resampling = T, resamplingSteps = 2, proposal = NULL, exponents = NULL, adaptive = T, proposalScale = 0.5, x=3.11, m=7E-08, sampling="multinomial",ess.limit=NULL,lastResample = 1,pars.lower=NULL, pars.upper=NULL,mutate="Metropolis",b=1e-04){ +smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, resampling = T, resamplingSteps = 2, proposal = NULL, exponents = NULL, adaptive = T, proposalScale = 0.5, x=3.11, m=7E-08, sampling="multinomial",ess.limit=NULL,lastResample = 1,pars.lower=NULL, pars.upper=NULL,mutate="Metropolis",b=1e-04, diagnostics = NULL){ if(resamplingSteps < 1) stop("SMC error, resamplingSteps can't be < 1") @@ -33,13 +34,10 @@ smcSampler <- function(bayesianSetup, reference=NULL, initialParticles = 1000, i info$dist.vec <- vector("numeric", iterations) info$res.of <- rep(FALSE, iterations) info$res.ess <- rep(FALSE, iterations) - - weights <- vector("numeric", initialParticles) - oldweights <- vector("numeric", initialParticles) - oldll <- vector("numeric", initialParticles) - oldInter <- vector("numeric", initialParticles) - expTrack <- vector("numeric", initialParticles) - + + + weights <- oldweights <- oldInter <- rep(0, initialParticles) + setup <- checkBayesianSetup(bayesianSetup) if(class(initialParticles) == "numeric"){ @@ -283,7 +281,7 @@ smcSampler <- function(bayesianSetup, reference=NULL, initialParticles = 1000, i } } } - info$dist.vec[(icount-1)] <- getSampleDistance(particles, getSample(reference),type="BH") + if(!is.null(diagnostics)) info$diagnostics[[icount-1]] <- diagnostics(particles) } info$rejectionRate = rejectionRate / (iterations * resamplingSteps) From f9f67331004ba89e90148cfe0e7b19b2011c73dc Mon Sep 17 00:00:00 2001 From: Florian Hartig Date: Thu, 9 Aug 2018 13:55:06 +0200 Subject: [PATCH 09/42] updates test --- Development/SMC/SMC-basicTests.R | 31 ++++++++++++++++++++++++++++--- 1 file changed, 28 insertions(+), 3 deletions(-) diff --git a/Development/SMC/SMC-basicTests.R b/Development/SMC/SMC-basicTests.R index 24d808f..cd9f065 100644 --- a/Development/SMC/SMC-basicTests.R +++ b/Development/SMC/SMC-basicTests.R @@ -1,14 +1,39 @@ -ll = function(x) sum(dnorm(x, log = T)) +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) - -out = smcSampler(bayesianSetup = setup, initialParticles = 2000, iterations = 1, ess.limit = 10000000) +out = smcSampler(bayesianSetup = setup, initialParticles = 2000, iterations = 2, ess.limit = 10000000, diagnostics = function(x)sum(mean(x))) hist(out$particles) 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 From 87fce94c43b94334f00a1dd33af02bae0706ad6c Mon Sep 17 00:00:00 2001 From: Florian Hartig Date: Fri, 10 Aug 2018 12:22:31 +0200 Subject: [PATCH 10/42] smc updates --- BayesianTools/R/SMC.R | 112 +++++++++++++++++------- BayesianTools/tests/testthat/test-SMC.R | 26 ++++++ Development/SMC/SMC-basicTests.R | 20 ++++- 3 files changed, 121 insertions(+), 37 deletions(-) create mode 100644 BayesianTools/tests/testthat/test-SMC.R diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index 9171740..32519ee 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -23,7 +23,27 @@ #' @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, exponents = NULL, adaptive = T, proposalScale = 0.5, x=3.11, m=7E-08, sampling="multinomial",ess.limit=NULL,lastResample = 1,pars.lower=NULL, pars.upper=NULL,mutate="Metropolis",b=1e-04, diagnostics = NULL){ +smcSampler <- function(bayesianSetup, + initialParticles = 1000, + iterations = 10, + resampling = T, + resamplingSteps = 2, + proposal = NULL, + exponents = NULL, + adaptive = T, + proposalScale = 0.5, + x=3.11, + m=7E-08, + sampling="multinomial", + ess.limit=NULL, + lastResample = 1, + pars.lower=NULL, + pars.upper=NULL, + mutate="Metropolis", + b=1e-04, + diagnostics = NULL){ + + ########### SETUP STEPS ######################## if(resamplingSteps < 1) stop("SMC error, resamplingSteps can't be < 1") @@ -36,12 +56,21 @@ smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, info$res.ess <- rep(FALSE, iterations) - weights <- oldweights <- oldInter <- rep(0, initialParticles) - setup <- checkBayesianSetup(bayesianSetup) + ### InitialParticles + #TODO documentation + if(class(initialParticles) == "numeric"){ initialParticles = bayesianSetup$prior$sampler(initialParticles) + importanceDensity = bayesianSetup$prior$density + } + if(class(initialParticles) == "matrix"){ + importanceDensity = bayesianSetup$prior$density + } + if(class(initialParticles) == "list"){ + initialParticles = initialParticles$particles + importanceDensity = initialParticles$density } if (any(is.infinite(setup$prior$density(initialParticles)))) stop("initialParticles outside prior range") @@ -51,6 +80,8 @@ smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, particleSize = nrow(initialParticles) if(is.null(ess.limit)){ess.limit <- round(particleSize) * 0.5} + weights <- oldweights <- oldInter <- rep(0, particleSize) + posterior = matrix(nrow = particleSize, ncol = 3) numPar <- ncol(initialParticles) @@ -59,7 +90,10 @@ smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, } usedUp = 0 + + ########### 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)){ @@ -90,6 +124,7 @@ smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, #for (i in 1:iterations){ while(icount <= length(exponents)){ + # Using a while loop instead of for because in the adaptive algorithm, the number of iterations may increase # Reweighting @@ -101,26 +136,22 @@ smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, # Jeremia is using beta sequence # https://github.com/florianhartig/BayesianTools/issues/23 - - posterior = setup$posterior$density(particles, returnAll = T) - priorValues <- posterior[,3] - likelihoodValues <- posterior[, 2] - postValues <- posterior[,1] - + print(c("weights", head(weights))) ## Update weights # Create intermediary distribution curExp <- exponents[icount] - interDist <- curExp * postValues + (1-curExp) * priorValues + interDist <- curExp * setup$posterior$density(particles) + + (1-curExp) * importanceDensity(particles) # Calculate new weights oldweights <- weights weights <- weights + (interDist - oldInter) # Normalize (log-)weights so that the sum (of non-logs) equals 1 - sumw <- sum(sort(exp(weights))) - weights <- weights - log(sumw) + + weights <- weights - BayesianTools:::logSumExp(weights) # It may occur that the difference between the new and old weights of a particle is so great that an overflow occurs. # In this case, the algorithm reduces the difference between the current and next intermediary distribution (i.e. the difference @@ -144,6 +175,7 @@ smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, # relativeL = exp(likelihoodValues - max(likelihoodValues, na.rm = T))^(1/iterations) + ######################################################### # Re?sampling step # Option 1) Always resample (default BT) @@ -167,31 +199,20 @@ smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, # on the location of particles in parameter space, the weights are not considered in the output) if(ess < ess.limit | icount == (length(exponents) - lastResample)){ oldExp <- curExp - if(sampling == "multinomial"){ - sel <- sample.int(n=length(likelihoodValues), size = length(likelihoodValues), replace = T, prob = exp(weights)) - } else if(sampling == "residual"){ - sel <- residualResampling(weights) - } else if(sampling == "systematic"){ - sel <- systematicResampling(weights) - } else{ - stop("Invalid string for resampling argument") - } - + + sel = resample(weights, method = sampling) + particles = particles[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 - sumw <- sum(sort(exp(weights))) - weights <- weights - log(sumw) + weights <- weights - BayesianTools:::logSumExp(weights) # Get posterior values and distribution for new particles - posterior = setup$posterior$density(particles, returnAll = T) - priorValues <- posterior[,3] - likelihoodValues <- posterior[, 2] - postValues <- posterior[,1] - oldInter <- oldExp * postValues + (1-oldExp) * priorValues + oldInter <- curExp * setup$posterior$density(particles) + + (1-curExp) * importanceDensity(particles) # Determine if exponents need to be modified if(ess < estar){ @@ -216,6 +237,7 @@ smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, if (numPar == 1) particles = matrix(particles, ncol = 1) + ########################################################### # Mutate # 1) When to do the resampling @@ -286,12 +308,13 @@ smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, info$rejectionRate = rejectionRate / (iterations * resamplingSteps) + settings = list(initialParticles = initialParticles, proposalGenerator = proposalGenerator) + out = list( setup = setup, - initialParticles = initialParticles, + settings = settings, particles = particles, - posteriorValues = posterior, - proposalGenerator = proposalGenerator, + posterior = posterior, info = info ) @@ -301,6 +324,9 @@ smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, } +#' Residual Resampling +#' +#' @keywords internal residualResampling <- function(weights){ weights <- exp(weights) # Define number of replications for each particle (integer), as well as residuals @@ -317,7 +343,6 @@ residualResampling <- function(weights){ new.parts <- rep(1:length(weights), nrep.int) new.parts <- c(new.parts, missing) return(new.parts) - } systematicResampling <- function(weights){ @@ -345,4 +370,23 @@ systematicResampling <- function(weights){ } return(new.parts) -} \ No newline at end of file +} + + + +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) +} + diff --git a/BayesianTools/tests/testthat/test-SMC.R b/BayesianTools/tests/testthat/test-SMC.R new file mode 100644 index 0000000..856e1b7 --- /dev/null +++ b/BayesianTools/tests/testthat/test-SMC.R @@ -0,0 +1,26 @@ +skip_on_cran() + +context("Test SMC and utility functions") + + + + + +test_that("resampling works", { + + weights = log(c(1,2,3)/6) + size = 10000 + + x1 = replicate(size, 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) + +} diff --git a/Development/SMC/SMC-basicTests.R b/Development/SMC/SMC-basicTests.R index cd9f065..7060437 100644 --- a/Development/SMC/SMC-basicTests.R +++ b/Development/SMC/SMC-basicTests.R @@ -1,6 +1,6 @@ ll = function(x){ - for(i in 1:10000) mean(rnorm(10)) + #for(i in 1:10000) mean(rnorm(10)) sum(dnorm(x, log = T)) } setup = createBayesianSetup(ll, lower = c(-10, -10), upper = c(10, 10)) @@ -10,8 +10,22 @@ out = runMCMC(bayesianSetup = setup, settings = settings) ks.test(getSample(out, start = 5000)[,1], msm::ptnorm, 0, 1 , -10, 10) plot(out) -out = smcSampler(bayesianSetup = setup, initialParticles = 2000, iterations = 2, ess.limit = 10000000, diagnostics = function(x)sum(mean(x))) -hist(out$particles) + +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) From d0a3d849ac61c2a832c4c75f8fdd73cd47ccf790 Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Thu, 16 Aug 2018 18:19:35 +0200 Subject: [PATCH 11/42] Fix some SMC issues - Recalculation of previous intermediary distribution for resampled particles - mutate as a separate function - extra mutate steps at the end of run --- BayesianTools/R/SMC.R | 377 +++++++++++++++++++++++++++++++----------- 1 file changed, 276 insertions(+), 101 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index 32519ee..5d42a15 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -28,6 +28,7 @@ smcSampler <- function(bayesianSetup, iterations = 10, resampling = T, resamplingSteps = 2, + lastMutateSteps = 5, # TODO document proposal = NULL, exponents = NULL, adaptive = T, @@ -39,8 +40,8 @@ smcSampler <- function(bayesianSetup, lastResample = 1, pars.lower=NULL, pars.upper=NULL, - mutate="Metropolis", - b=1e-04, + mutate.method ="Metropolis", # TODO document + b=1e-04, # TODO document diagnostics = NULL){ ########### SETUP STEPS ######################## @@ -55,7 +56,6 @@ smcSampler <- function(bayesianSetup, info$res.of <- rep(FALSE, iterations) info$res.ess <- rep(FALSE, iterations) - setup <- checkBayesianSetup(bayesianSetup) ### InitialParticles @@ -69,10 +69,12 @@ smcSampler <- function(bayesianSetup, importanceDensity = bayesianSetup$prior$density } if(class(initialParticles) == "list"){ - initialParticles = initialParticles$particles importanceDensity = initialParticles$density + initialParticles = as.matrix(initialParticles$particles,ncol=1) } + #print(c("density", setup$prior$density(initialParticles))) + if (any(is.infinite(setup$prior$density(initialParticles)))) stop("initialParticles outside prior range") particles <- initialParticles @@ -85,7 +87,7 @@ smcSampler <- function(bayesianSetup, posterior = matrix(nrow = particleSize, ncol = 3) numPar <- ncol(initialParticles) - if(mutate %in% c("Metropolis", "adaptive")){ + if(mutate.method %in% c("Metropolis", "adaptive")){ if (is.null(proposal)) proposalGenerator = createProposalGenerator(rep(40,numPar)) } @@ -102,6 +104,8 @@ smcSampler <- function(bayesianSetup, exponents <- pmin(1,exponents) exponents <- exponents[2:length(exponents)] } + print(c("length xponents", length(exponents))) + print(c("exponents", 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 @@ -114,19 +118,128 @@ smcSampler <- function(bayesianSetup, # All particles are given equal weight at initialization weights[1:length(weights)] <- log(1/particleSize) + oldweights <- weights + + # Initial importance distribution + importanceValues <- importanceDensity(particles) + posteriorValues <- setup$posterior$density(particles) + #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, message = TRUE)} + if (is.null(proposal)){ proposalGenerator = createProposalGenerator(sds)} + + + ########################################################## + # Loop starts here # + ########################################################## #for (i in 1:iterations){ while(icount <= length(exponents)){ + #posteriorValues <- setup$posterior$density(particles) + #importanceValues <- importanceDensity(particles) + + if (numPar == 1) particles = matrix(particles, ncol = 1) + + print(c("icount",icount)) + print(c("posteriorValues",head(posteriorValues))) + print(c("importanceValues",head(importanceValues))) + # 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 + + #if (resampling == T){ + + + ####------------- Replaced by mutate function - Start ------------------- + # if(mutate=="Metropolis"){ + # + # if (adaptive == T){ + # proposalGenerator = updateProposalGenerator(proposalGenerator, particles) + # } + # + # for (j in 1:resamplingSteps){ + # particlesProposals = proposalGenerator$returnProposalMatrix(particles, scale = proposalScale) + # proposalPosteriors <- setup$posterior$density(particlesProposals) + # proposalImportance <- importanceDensity(particlesProposals) + # + # #jumpProb <- exp(setup$posterior$density(particlesProposals) - likelihoodValues[sel])^(i/iterations) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) + # #jumpProb <- exp(setup$posterior$density(particlesProposals) - setup$posterior$density(particles)) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) + # jumpProb <- exp(proposalPosteriors - posteriorValues) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) + # #jumpProb <- exp(setup$posterior$density(particlesProposals) - posteriorValues) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) + # + # + # print(c("particlesProposals", head(particlesProposals))) + # print(c("particles", head(particles))) + # accepted <- jumpProb > runif(length(jumpProb), 0 ,1) + # rejectionRate = rejectionRate + sum(accepted) + # particles[accepted, ] = particlesProposals[accepted, ] + # posteriorValues[accepted] <- proposalPosteriors[accepted] + # importanceValues[accepted] <- proposalImportance[accepted] + # info$resamplingAcceptance[(icount),j] <- sum(accepted)/particleSize + # } + # } else if(mutate=="DE"){ + # # Differential evolution + # # For now (test): basic DE algorithm + # + # for(j in 1:resamplingSteps){ + # + # # Loop over particles + # for(part in 1:particleSize){ + # particlesOld <- particles + # particleDiff <- rep(0,numPar) + # newParts <- rep(0L,2) + # + # # Sample 2 other particles + # while(all(particleDiff==0) | part %in% newParts){ + # # 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=particleSize,size=2) + # particleDiff <- particlesOld[newParts[1],] - particlesOld[newParts[2],] + # } + # + # particleProp <- particlesOld[part,] + (particleDiff * proposalScale) + runif(numPar,-b,b) + # jumpProb <- exp(setup$posterior$density(particleProp) - likelihoodValues[part]) * exp(setup$prior$density(particleProp) - setup$prior$density(particles[part,])) + # + # if(jumpProb > runif(length(jumpProb), 0 ,1)){ + # particles[part,] <- particleProp + # if(is.na(info$resamplingAcceptance[(icount-1),j])){ + # info$resamplingAcceptance[(icount-1),j] <- 1/particleSize + # } else{ + # info$resamplingAcceptance[(icount-1),j] <- info$resamplingAcceptance[(icount-1),j] + 1/particleSize + # } + # } else{ + # particles[part,] <- particlesOld[part,] + # } + # } + # } + # } + ####------------- Replaced by mutate function - End ------------------- + + mutate.out <- mutate(setup = setup, particles = particles, proposalGenerator = proposalGenerator, posteriorValues = posteriorValues, importanceDensity = importanceDensity, method = mutate.method, steps = resamplingSteps, proposalScale = proposalScale, adaptive = adaptive) + 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 @@ -136,17 +249,16 @@ smcSampler <- function(bayesianSetup, # Jeremia is using beta sequence # https://github.com/florianhartig/BayesianTools/issues/23 - - print(c("weights", head(weights))) ## Update weights # Create intermediary distribution curExp <- exponents[icount] - interDist <- curExp * setup$posterior$density(particles) - + (1-curExp) * importanceDensity(particles) + + oldInter <- oldExp * posteriorValues + (1-oldExp) * importanceValues + interDist <- curExp * posteriorValues + (1-curExp) * importanceValues # Calculate new weights - oldweights <- weights + weights <- weights + (interDist - oldInter) # Normalize (log-)weights so that the sum (of non-logs) equals 1 @@ -168,6 +280,16 @@ smcSampler <- function(bayesianSetup, next } + print("-------") + print(c("icount", icount)) + print(c("particles", head(particles))) + print(c("importance", head(importanceValues))) + print(c("posterior", head(posteriorValues))) + print(c("interDist", head(interDist))) + print(c("oldInter", head(oldInter))) + print(c("weights", head(weights))) + print(c("oldweights", head(oldweights))) + # posterior = setup$posterior$density(particles, returnAll = T) # likelihoodValues <- posterior[,2] # ll = likelihoodValues - max(likelihoodValues, na.rm = T) @@ -188,31 +310,19 @@ smcSampler <- function(bayesianSetup, # Calculate effective sample size (ESS) ess <- 1 / sum(exp(2 * weights)) - - # Save intermediate distribution for the next iteration - oldInter <- interDist if(is.infinite(ess)){print(weights)} # 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 | icount == (length(exponents) - lastResample)){ - oldExp <- curExp + + print(c("ess", ess)) + print(c("ess.limit", ess.limit)) + print(c("curExp", curExp)) - sel = resample(weights, method = sampling) - - particles = particles[sel,] - info$survivingParticles[icount] = length(unique(sel)) + if(ess < ess.limit | icount == (length(exponents) - lastResample)){ - # 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) - - # Get posterior values and distribution for new particles - oldInter <- curExp * setup$posterior$density(particles) - + (1-curExp) * importanceDensity(particles) # Determine if exponents need to be modified if(ess < estar){ @@ -223,89 +333,75 @@ smcSampler <- function(bayesianSetup, expdiff <- (exponents[icount] - exponents[(icount-1)]) * 0.5 newExp <- exponents[(icount-1)] + expdiff - exponents <- append(exponents,newExp,after=(icount-1)) + print(c("exponents[icount]",exponents[icount])) + print(c("exponents[icount-1]",exponents[icount-1])) + print(c("expdiff", expdiff)) + print(c("newExp",newExp)) + + weights <- oldweights + + if(icount==1){ + newExp <- exponents[icount] * 0.5 + exponents <- c(newExp,exponents) + } else { + exponents <- append(exponents,newExp,after=(icount-1)) + } + + print(c("length(exponents", length(exponents))) + print(c("new exponents", exponents)) info$res.ess[icount] <- TRUE next } + + oldExp <- curExp + + sel = resample(weights, method = sampling) + + particles = particles[sel,] + 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) + + # Get posterior values and distribution for new particles + #oldInter <- curExp * setup$posterior$density(particles) + # + (1-curExp) * importanceDensity(particles) + #oldInter <- curExp * posteriorValues + (1-curExp) * importanceValues + + print("-- Resampling") + print(c("posterior", head(posteriorValues))) + print(c("importance", head(importanceValues))) + #print(c("oldInter", head(oldInter))) + + + } else{ + # Save intermediate distribution for the next iteration + #oldInter <- interDist } + oldweights <- weights + + #print(c("oldInter saved", head(oldInter))) + info$ess.vec[icount] <- ess icount <- icount + 1 - curExp <- m * (icount*200/iterations)^x + #curExp <- m * (icount*200/iterations)^x oldExp <- curExp - if (numPar == 1) particles = matrix(particles, ncol = 1) - - ########################################################### - # 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 - - #if (resampling == T){ - if(mutate=="Metropolis"){ - - if (adaptive == T){ - proposalGenerator = updateProposalGenerator(proposalGenerator, particles) - } - - 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)) - jumpProb <- exp(setup$posterior$density(particlesProposals) - likelihoodValues) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) - - accepted <- jumpProb > runif(length(jumpProb), 0 ,1) - rejectionRate = rejectionRate + sum(accepted) - particles[accepted, ] = particlesProposals[accepted, ] - info$resamplingAcceptance[(icount-1),j] <- sum(accepted)/particleSize - } - } else if(mutate=="DE"){ - # Differential evolution - # For now (test): basic DE algorithm - - for(j in 1:resamplingSteps){ - - # Loop over particles - for(part in 1:particleSize){ - particlesOld <- particles - particleDiff <- rep(0,numPar) - newParts <- rep(0L,2) - - # Sample 2 other particles - while(all(particleDiff==0) | part %in% newParts){ - # 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=particleSize,size=2) - particleDiff <- particlesOld[newParts[1],] - particlesOld[newParts[2],] - } - - particleProp <- particlesOld[part,] + (particleDiff * proposalScale) + runif(numPar,-b,b) - jumpProb <- exp(setup$posterior$density(particleProp) - likelihoodValues[part]) * exp(setup$prior$density(particleProp) - setup$prior$density(particles[part,])) - - if(jumpProb > runif(length(jumpProb), 0 ,1)){ - particles[part,] <- particleProp - if(is.na(info$resamplingAcceptance[(icount-1),j])){ - info$resamplingAcceptance[(icount-1),j] <- 1/particleSize - } else{ - info$resamplingAcceptance[(icount-1),j] <- info$resamplingAcceptance[(icount-1),j] + 1/particleSize - } - } else{ - particles[part,] <- particlesOld[part,] - } - } - } - } if(!is.null(diagnostics)) info$diagnostics[[icount-1]] <- diagnostics(particles) } + # Last mutation, to increase diversity between particles + mutate.out <- mutate(setup = setup, particles = particles, proposalGenerator = proposalGenerator, posteriorValues = posteriorValues, importanceDensity = importanceDensity, method = mutate.method, steps = lastMutateSteps, proposalScale = proposalScale, adaptive = adaptive) + particles <- mutate.out$particles + posteriorValues <- mutate.out$posteriorValues + importanceValues <- mutate.out$importanceValues + info$rejectionRate = rejectionRate / (iterations * resamplingSteps) settings = list(initialParticles = initialParticles, proposalGenerator = proposalGenerator) @@ -352,7 +448,7 @@ systematicResampling <- function(weights){ rank.weights <- rank(weights, ties.method = "first") sort.weights <- sort(weights) cumu.weights <- c(0,cumsum(sort.weights)[1:(length(weights)-1)]) - print(cumu.weights) + #print(cumu.weights) n.parts <- length(weights) # Number of particles u <- runif(n=1, min=0, max=1/n.parts) @@ -361,10 +457,10 @@ systematicResampling <- function(weights){ new.parts.sort <- vector("numeric", n.parts) for(i in 1:n.parts){ U[i] <- ((i-1)/n.parts) + u - print(c(i, U[i])) + #print(c(i, U[i])) new.parts.sort[i] <- tail(which(cumu.weights <= U[i]),1) - print(c("new.parts.sort", new.parts.sort[i])) + #print(c("new.parts.sort", new.parts.sort[i])) # Match sorted ranks back to original order of weights new.parts[i] <- which(rank.weights==new.parts.sort[i]) } @@ -390,3 +486,82 @@ resample <- function(weights, method = "multinomial"){ return(sel) } +beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldInter, curWeights, curExp, tol=0.001){ + # A function to dynamically set the next exponent to build the next intermediary distribution. + # Uses the bisection method. Following Jasra et al. (2001), Scand J Statist, doi: 10.1111/j.1467-9469.2010.00723.x + + # Initial exponent - set to 1 (maximum possible value) + tryDiff <- 1 - curExp + + while(abs(ess-target.ess) > tol){ + tryExp <- curExp + tryDiff + + tryDist <- tryExp * posteriorValues + (1-tryExp) * importanceValues + + tryWeights <- curWeights + (tryDist - oldInter) + # Normalize (log-)weights so that the sum (of non-logs) equals 1 + tryWeights <- tryWeights - BayesianTools:::logSumExp(tryWeights) + + try.ess <- 1 / sum(exp(2 * weights)) + + if(try.ess - target.ess > tol){ + # Greater ESS than desired -> choose a *larger* exponent in the next iteration (efficiency) + tryDiff <- tryDiff + tryDiff * 0.5 + + # Exponent may not be greater than one + if(tryDiff + curExp >= 1){ + tryExp <- 1 + break + } + + } 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) + tryDiff <- tryDiff - (tryDiff * 0.5) + + # Exponent may not decrease from previous value + if(tryDiff < 0){ + # This should not occur - leaving it for testing only + stop("problem defining next exponent") + } + } + } + + out <- list(newExp=tryExp, weights=tryWeights, interDist=tryDist, ess=try.ess) +} + + +mutate <- function(setup, particles, proposalGenerator, posteriorValues, importanceDensity, method, steps, proposalScale, adaptive = TRUE){ + if(is.vector(particles)){particles = matrix(particles, ncol = 1)} + acceptance <- vector("numeric", length=steps) + importanceValues <- importanceDensity(particles) + + if(method=="Metropolis"){ + if(adaptive){ + proposalGenerator = updateProposalGenerator(proposalGenerator, particles) + } + + for(j in 1:steps){ + particlesProposals = proposalGenerator$returnProposalMatrix(particles, scale = proposalScale) + proposalPosteriors <- setup$posterior$density(particlesProposals) + proposalImportance <- importanceDensity(particlesProposals) + + #jumpProb <- exp(setup$posterior$density(particlesProposals) - likelihoodValues[sel])^(i/iterations) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) + #jumpProb <- exp(setup$posterior$density(particlesProposals) - setup$posterior$density(particles)) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) + jumpProb <- exp(proposalPosteriors - posteriorValues) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) + #jumpProb <- exp(setup$posterior$density(particlesProposals) - posteriorValues) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) + + + print(c("particlesProposals", head(particlesProposals))) + print(c("particles", head(particles))) + accepted <- jumpProb > runif(length(jumpProb), 0 ,1) + #rejectionRate = rejectionRate + sum(accepted) + particles[accepted, ] = particlesProposals[accepted, ] + posteriorValues[accepted] <- proposalPosteriors[accepted] + importanceValues[accepted] <- proposalImportance[accepted] + acceptance[j] <- sum(accepted)/nrow(particles) + #info$resamplingAcceptance[(icount),j] <- sum(accepted)/particleSize + } + } + out <- list(particles=particles, posteriorValues = posteriorValues, importanceValues = importanceValues, acceptance = acceptance) +} From 96086e4691353c3f49b65a3d3ee27201b940bf52 Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Mon, 20 Aug 2018 13:45:18 +0200 Subject: [PATCH 12/42] Fix algorithm for adaptive sequence --- BayesianTools/R/SMC.R | 192 +++++++++++++++++++++++++----------------- 1 file changed, 115 insertions(+), 77 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index 5d42a15..fb40678 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -37,6 +37,7 @@ smcSampler <- function(bayesianSetup, m=7E-08, sampling="multinomial", ess.limit=NULL, + ess.factor = 0.95, # TODO document lastResample = 1, pars.lower=NULL, pars.upper=NULL, @@ -98,21 +99,20 @@ smcSampler <- function(bayesianSetup, # 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)] - } - print(c("length xponents", length(exponents))) - print(c("exponents", exponents)) + # 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 <- exponents[1] + curExp <- 0 # For fully adaptive algorithm: initial value of curExp does not matter, but must be < 1 icount <- 1 # Iterations counter @@ -133,22 +133,26 @@ smcSampler <- function(bayesianSetup, } if (is.null(proposal)){ proposalGenerator = createProposalGenerator(sds)} + lastIteration <- FALSE + ########################################################## # Loop starts here # ########################################################## #for (i in 1:iterations){ - while(icount <= length(exponents)){ + # while(icount <= length(exponents)){ + while(curExp <= 1 & !lastIteration){ # For fully adaptive algorithm + if(curExp == 1){lastIteration <- TRUE} #posteriorValues <- setup$posterior$density(particles) #importanceValues <- importanceDensity(particles) if (numPar == 1) particles = matrix(particles, ncol = 1) - print(c("icount",icount)) - print(c("posteriorValues",head(posteriorValues))) - print(c("importanceValues",head(importanceValues))) + # print(c("icount",icount)) + # print(c("posteriorValues",head(posteriorValues))) + # print(c("importanceValues",head(importanceValues))) # Using a while loop instead of for because in the adaptive algorithm, the number of iterations may increase @@ -250,45 +254,54 @@ smcSampler <- function(bayesianSetup, # https://github.com/florianhartig/BayesianTools/issues/23 + ## Update weights # Create intermediary distribution - curExp <- exponents[icount] - + #curExp <- exponents[icount] + ess <- 1 / sum(exp(2 * weights)) oldInter <- oldExp * posteriorValues + (1-oldExp) * importanceValues - interDist <- curExp * posteriorValues + (1-curExp) * importanceValues + inter.out <- beta.search(ess=ess, target.ess = (ess * ess.factor), 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 + + #interDist <- curExp * posteriorValues + (1-curExp) * importanceValues # Calculate new weights weights <- weights + (interDist - oldInter) - # Normalize (log-)weights so that the sum (of non-logs) equals 1 - weights <- weights - BayesianTools:::logSumExp(weights) + + + # It may occur that the difference between the new and old weights of a particle is so great that an overflow occurs. # In this case, the algorithm reduces the difference between the current and next intermediary distribution (i.e. the difference # between the current and next value of the exponent). As a result, one more iteration is added. # This is the same that happens when the effective sample size becomes *very* low. #if(is.infinite(sum(sort(exp(weights))))){ - if(any(is.infinite(weights))){ - print("Infinite weights sum; retrying with new intermediary distribution") - expdiff <- (exponents[icount] - exponents[(icount-1)]) * 0.5 - newExp <- exponents[(icount-1)] + expdiff - exponents <- append(exponents,newExp,after=(icount-1)) - weights <- oldweights - info$res.of[icount] <- TRUE - next - } + + # if(any(is.infinite(weights))){ + # print("Infinite weights sum; retrying with new intermediary distribution") + # expdiff <- (exponents[icount] - exponents[(icount-1)]) * 0.5 + # newExp <- exponents[(icount-1)] + expdiff + # exponents <- append(exponents,newExp,after=(icount-1)) + # weights <- oldweights + # info$res.of[icount] <- TRUE + # next + # } print("-------") print(c("icount", icount)) - print(c("particles", head(particles))) - print(c("importance", head(importanceValues))) - print(c("posterior", head(posteriorValues))) - print(c("interDist", head(interDist))) - print(c("oldInter", head(oldInter))) - print(c("weights", head(weights))) - print(c("oldweights", head(oldweights))) + # print(c("particles", head(particles))) + # print(c("importance", head(importanceValues))) + # print(c("posterior", head(posteriorValues))) + # print(c("interDist", head(interDist))) + # print(c("oldInter", head(oldInter))) + # print(c("weights", head(weights))) + # print(c("oldweights", head(oldweights))) # posterior = setup$posterior$density(particles, returnAll = T) # likelihoodValues <- posterior[,2] @@ -309,7 +322,7 @@ smcSampler <- function(bayesianSetup, #sel = sample.int(n=length(likelihoodValues), size = length(likelihoodValues), replace = T, prob = relativeL) # Calculate effective sample size (ESS) - ess <- 1 / sum(exp(2 * weights)) + #ess <- 1 / sum(exp(2 * weights)) if(is.infinite(ess)){print(weights)} @@ -325,33 +338,33 @@ smcSampler <- function(bayesianSetup, # Determine if exponents need to be modified - if(ess < estar){ - # If yes: - # - discard newly sampled population of particles - # - do not increase counter - # - modify x - - expdiff <- (exponents[icount] - exponents[(icount-1)]) * 0.5 - newExp <- exponents[(icount-1)] + expdiff - print(c("exponents[icount]",exponents[icount])) - print(c("exponents[icount-1]",exponents[icount-1])) - print(c("expdiff", expdiff)) - print(c("newExp",newExp)) - - weights <- oldweights - - if(icount==1){ - newExp <- exponents[icount] * 0.5 - exponents <- c(newExp,exponents) - } else { - exponents <- append(exponents,newExp,after=(icount-1)) - } - - print(c("length(exponents", length(exponents))) - print(c("new exponents", exponents)) - info$res.ess[icount] <- TRUE - next - } + # if(ess < estar){ + # # If yes: + # # - discard newly sampled population of particles + # # - do not increase counter + # # - modify x + # + # expdiff <- (exponents[icount] - exponents[(icount-1)]) * 0.5 + # newExp <- exponents[(icount-1)] + expdiff + # print(c("exponents[icount]",exponents[icount])) + # print(c("exponents[icount-1]",exponents[icount-1])) + # print(c("expdiff", expdiff)) + # print(c("newExp",newExp)) + # + # weights <- oldweights + # + # if(icount==1){ + # newExp <- exponents[icount] * 0.5 + # exponents <- c(newExp,exponents) + # } else { + # exponents <- append(exponents,newExp,after=(icount-1)) + # } + # + # print(c("length(exponents", length(exponents))) + # print(c("new exponents", exponents)) + # info$res.ess[icount] <- TRUE + # next + # } oldExp <- curExp @@ -373,8 +386,8 @@ smcSampler <- function(bayesianSetup, #oldInter <- curExp * posteriorValues + (1-curExp) * importanceValues print("-- Resampling") - print(c("posterior", head(posteriorValues))) - print(c("importance", head(importanceValues))) + # print(c("posterior", head(posteriorValues))) + # print(c("importance", head(importanceValues))) #print(c("oldInter", head(oldInter))) @@ -419,6 +432,8 @@ smcSampler <- function(bayesianSetup, } +###################################################### +# Auxiliary functions for resampling #' Residual Resampling #' @@ -486,15 +501,24 @@ resample <- function(weights, method = "multinomial"){ return(sel) } -beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldInter, curWeights, curExp, tol=0.001){ +###################################################### + +# Auxiliary function for adaptive algorithm + +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. (2001), 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(ess-target.ess) > tol){ - tryExp <- curExp + tryDiff + while(abs(try.ess-target.ess) > tol & tryExp < 1){ + #tryExp <- curExp + tryDiff + tryExp <- (a+b) * 0.5 tryDist <- tryExp * posteriorValues + (1-tryExp) * importanceValues @@ -502,34 +526,48 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI # Normalize (log-)weights so that the sum (of non-logs) equals 1 tryWeights <- tryWeights - BayesianTools:::logSumExp(tryWeights) - try.ess <- 1 / sum(exp(2 * weights)) + try.ess <- 1 / sum(exp(2 * tryWeights)) + + + print("---------------") + print(c("ess", ess)) + print(c("target.ess", target.ess)) + print(c("try.ess", try.ess)) + print(c("a,b", a,b)) if(try.ess - target.ess > tol){ # Greater ESS than desired -> choose a *larger* exponent in the next iteration (efficiency) - tryDiff <- tryDiff + tryDiff * 0.5 + #tryDiff <- tryDiff + tryDiff * 0.5 + a <- tryExp # Exponent may not be greater than one - if(tryDiff + curExp >= 1){ - tryExp <- 1 - break - } + #if(tryDiff + curExp >= 1){ + # tryExp <- 1 + # break + #} } 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) - tryDiff <- tryDiff - (tryDiff * 0.5) + #tryDiff <- tryDiff - (tryDiff * 0.5) + + b <- tryExp # Exponent may not decrease from previous value - if(tryDiff < 0){ + #if(tryDiff < 0){ # This should not occur - leaving it for testing only - stop("problem defining next exponent") - } + # stop("problem defining next exponent") + #} } } out <- list(newExp=tryExp, weights=tryWeights, interDist=tryDist, ess=try.ess) + return(out) } +###################################################### +# Auxiliary function for particle mutation + mutate <- function(setup, particles, proposalGenerator, posteriorValues, importanceDensity, method, steps, proposalScale, adaptive = TRUE){ if(is.vector(particles)){particles = matrix(particles, ncol = 1)} From 93becbee762650fc890a117a31f6df530cfbae17 Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Tue, 28 Aug 2018 10:28:00 +0200 Subject: [PATCH 13/42] Fix differential evolution mutate step --- BayesianTools/R/SMC.R | 251 +++++++++++++++++++----------------------- 1 file changed, 115 insertions(+), 136 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index fb40678..34e3477 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -50,12 +50,12 @@ smcSampler <- function(bayesianSetup, if(resamplingSteps < 1) stop("SMC error, resamplingSteps can't be < 1") info = list() - info$resamplingAcceptance = as.data.frame(matrix(nrow = iterations, ncol = resamplingSteps)) #Using DF instead of matrix because number of iterations may change due to adaptive algorithm - info$survivingParticles = rep(NA, iterations) - info$ess.vec <- vector("numeric", iterations) - info$dist.vec <- vector("numeric", iterations) - info$res.of <- rep(FALSE, iterations) - info$res.ess <- rep(FALSE, iterations) + # 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). + 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) + setup <- checkBayesianSetup(bayesianSetup) @@ -74,8 +74,6 @@ smcSampler <- function(bayesianSetup, initialParticles = as.matrix(initialParticles$particles,ncol=1) } - #print(c("density", setup$prior$density(initialParticles))) - if (any(is.infinite(setup$prior$density(initialParticles)))) stop("initialParticles outside prior range") particles <- initialParticles @@ -142,17 +140,12 @@ smcSampler <- function(bayesianSetup, #for (i in 1:iterations){ # while(icount <= length(exponents)){ - while(curExp <= 1 & !lastIteration){ # For fully adaptive algorithm - if(curExp == 1){lastIteration <- TRUE} - - #posteriorValues <- setup$posterior$density(particles) - #importanceValues <- importanceDensity(particles) + #while(curExp <= 1 & !lastIteration){ # For fully adaptive algorithm + # if(curExp == 1){lastIteration <- TRUE} + while(curExp < 1){ if (numPar == 1) particles = matrix(particles, ncol = 1) - # print(c("icount",icount)) - # print(c("posteriorValues",head(posteriorValues))) - # print(c("importanceValues",head(importanceValues))) # Using a while loop instead of for because in the adaptive algorithm, the number of iterations may increase @@ -236,7 +229,7 @@ smcSampler <- function(bayesianSetup, # } ####------------- Replaced by mutate function - End ------------------- - mutate.out <- mutate(setup = setup, particles = particles, proposalGenerator = proposalGenerator, posteriorValues = posteriorValues, importanceDensity = importanceDensity, method = mutate.method, steps = resamplingSteps, proposalScale = proposalScale, adaptive = adaptive) + 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 @@ -257,7 +250,6 @@ smcSampler <- function(bayesianSetup, ## Update weights # Create intermediary distribution - #curExp <- exponents[icount] ess <- 1 / sum(exp(2 * weights)) oldInter <- oldExp * posteriorValues + (1-oldExp) * importanceValues inter.out <- beta.search(ess=ess, target.ess = (ess * ess.factor), posteriorValues = posteriorValues, importanceValues = importanceValues, oldInter = oldInter, curWeights = weights, curExp = curExp) @@ -265,16 +257,22 @@ smcSampler <- function(bayesianSetup, 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 + } - #interDist <- curExp * posteriorValues + (1-curExp) * importanceValues + print(c("weights U", head(weights))) + print(c("weights sum", BayesianTools:::logSumExp(weights))) + print(c("interDist", head(interDist))) + print(c("oldInter", head(oldInter))) # Calculate new weights - - weights <- weights + (interDist - oldInter) + #weights <- weights + (interDist - oldInter) # Normalize (log-)weights so that the sum (of non-logs) equals 1 - weights <- weights - BayesianTools:::logSumExp(weights) - - + #weights <- weights - BayesianTools:::logSumExp(weights) # It may occur that the difference between the new and old weights of a particle is so great that an overflow occurs. @@ -295,19 +293,6 @@ smcSampler <- function(bayesianSetup, print("-------") print(c("icount", icount)) - # print(c("particles", head(particles))) - # print(c("importance", head(importanceValues))) - # print(c("posterior", head(posteriorValues))) - # print(c("interDist", head(interDist))) - # print(c("oldInter", head(oldInter))) - # print(c("weights", head(weights))) - # print(c("oldweights", head(oldweights))) - - # posterior = setup$posterior$density(particles, returnAll = T) - # likelihoodValues <- posterior[,2] - # ll = likelihoodValues - max(likelihoodValues, na.rm = T) - # llCutoff = sort(ll)[acceptanceTarget] - # relativeL = exp(likelihoodValues - max(likelihoodValues, na.rm = T))^(1/iterations) ######################################################### @@ -318,11 +303,7 @@ smcSampler <- function(bayesianSetup, # => for now, option 1 is activated if the argument ess.limit is set to a value >= the number of particles # Question: last step complete resampling, so that all particles are the same? - - #sel = sample.int(n=length(likelihoodValues), size = length(likelihoodValues), replace = T, prob = relativeL) - - # Calculate effective sample size (ESS) - #ess <- 1 / sum(exp(2 * weights)) + if(is.infinite(ess)){print(weights)} @@ -334,38 +315,8 @@ smcSampler <- function(bayesianSetup, print(c("ess.limit", ess.limit)) print(c("curExp", curExp)) - if(ess < ess.limit | icount == (length(exponents) - lastResample)){ - - - # Determine if exponents need to be modified - # if(ess < estar){ - # # If yes: - # # - discard newly sampled population of particles - # # - do not increase counter - # # - modify x - # - # expdiff <- (exponents[icount] - exponents[(icount-1)]) * 0.5 - # newExp <- exponents[(icount-1)] + expdiff - # print(c("exponents[icount]",exponents[icount])) - # print(c("exponents[icount-1]",exponents[icount-1])) - # print(c("expdiff", expdiff)) - # print(c("newExp",newExp)) - # - # weights <- oldweights - # - # if(icount==1){ - # newExp <- exponents[icount] * 0.5 - # exponents <- c(newExp,exponents) - # } else { - # exponents <- append(exponents,newExp,after=(icount-1)) - # } - # - # print(c("length(exponents", length(exponents))) - # print(c("new exponents", exponents)) - # info$res.ess[icount] <- TRUE - # next - # } - + if(ess < ess.limit | icount == (length(exponents) - lastResample) | doResample){ + oldExp <- curExp sel = resample(weights, method = sampling) @@ -380,43 +331,42 @@ smcSampler <- function(bayesianSetup, # Normalize (log-)weights so that the sum (of non-logs) equals 1 weights <- weights - BayesianTools:::logSumExp(weights) - # Get posterior values and distribution for new particles - #oldInter <- curExp * setup$posterior$density(particles) - # + (1-curExp) * importanceDensity(particles) - #oldInter <- curExp * posteriorValues + (1-curExp) * importanceValues - print("-- Resampling") - # print(c("posterior", head(posteriorValues))) - # print(c("importance", head(importanceValues))) - #print(c("oldInter", head(oldInter))) - - - } else{ - # Save intermediate distribution for the next iteration - #oldInter <- interDist - } + 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 - - #print(c("oldInter saved", head(oldInter))) - info$ess.vec[icount] <- ess icount <- icount + 1 - #curExp <- m * (icount*200/iterations)^x oldExp <- curExp if(!is.null(diagnostics)) info$diagnostics[[icount-1]] <- diagnostics(particles) } + # Last resampling step, so that particles are distributed according to target distribution + sel = resample(weights, method = sampling) + particles = particles[sel,] + # Last mutation, to increase diversity between particles - mutate.out <- mutate(setup = setup, particles = particles, proposalGenerator = proposalGenerator, posteriorValues = posteriorValues, importanceDensity = importanceDensity, method = mutate.method, steps = lastMutateSteps, proposalScale = proposalScale, adaptive = adaptive) + mutate.out <- mutate(setup = setup, particles = particles, proposalGenerator = proposalGenerator, posteriorValues = posteriorValues, importanceDensity = importanceDensity, method = mutate.method, steps = lastMutateSteps, proposalScale = proposalScale, adaptive = adaptive, b=b) particles <- mutate.out$particles posteriorValues <- mutate.out$posteriorValues importanceValues <- mutate.out$importanceValues + lastAccept <- mutate.out$acceptance info$rejectionRate = rejectionRate / (iterations * resamplingSteps) + info$resamplingAcceptance <- info$resamplingAcceptance[1:(icount-1),] + info$survivingParticles <- info$survivingParticles[1:(icount-1)] + info$exponents <- info$exponents[1:(icount-1)] + info$ess.vec <- info$ess.vec[1:(icount-1)] + info$lastAccept <- lastAccept + settings = list(initialParticles = initialParticles, proposalGenerator = proposalGenerator) out = list( @@ -457,13 +407,11 @@ residualResampling <- function(weights){ } systematicResampling <- function(weights){ - #print(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)]) - #print(cumu.weights) n.parts <- length(weights) # Number of particles u <- runif(n=1, min=0, max=1/n.parts) @@ -472,10 +420,9 @@ systematicResampling <- function(weights){ new.parts.sort <- vector("numeric", n.parts) for(i in 1:n.parts){ U[i] <- ((i-1)/n.parts) + u - #print(c(i, U[i])) + new.parts.sort[i] <- tail(which(cumu.weights <= U[i]),1) - #print(c("new.parts.sort", new.parts.sort[i])) # Match sorted ranks back to original order of weights new.parts[i] <- which(rank.weights==new.parts.sort[i]) } @@ -516,10 +463,20 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI 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){ - #tryExp <- curExp + tryDiff + print("*-*-*_*-*-*") + print(c("ess", ess)) + print(c("curWeights", head(curWeights))) + + while(abs(try.ess-target.ess) > tol & tryExp < 1 & abs(b-a) > 1e-10){ tryExp <- (a+b) * 0.5 + # print(c("ess", ess)) + # print(c("target.ess", target.ess)) + # print(c("try.ess", try.ess)) + # print(c("curExp", curExp)) + # print(c("tryExp", tryExp)) + # print(c("a,b", a, b)) + tryDist <- tryExp * posteriorValues + (1-tryExp) * importanceValues tryWeights <- curWeights + (tryDist - oldInter) @@ -528,40 +485,30 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI try.ess <- 1 / sum(exp(2 * tryWeights)) + #print(c("new try.ess", try.ess)) - print("---------------") - print(c("ess", ess)) - print(c("target.ess", target.ess)) - print(c("try.ess", try.ess)) - print(c("a,b", a,b)) - - if(try.ess - target.ess > tol){ + if(try.ess - target.ess > tol & !any(is.infinite(tryWeights))){ # Greater ESS than desired -> choose a *larger* exponent in the next iteration (efficiency) - #tryDiff <- tryDiff + tryDiff * 0.5 a <- tryExp - - # Exponent may not be greater than one - #if(tryDiff + curExp >= 1){ - # tryExp <- 1 - # break - #} - } 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) - #tryDiff <- tryDiff - (tryDiff * 0.5) - b <- tryExp - - # Exponent may not decrease from previous value - #if(tryDiff < 0){ - # This should not occur - leaving it for testing only - # stop("problem defining next exponent") - #} } } - out <- list(newExp=tryExp, weights=tryWeights, interDist=tryDist, ess=try.ess) + if(any(is.infinite(tryWeights)) | is.infinite(try.ess)){ + doResample <- TRUE + tryExp <- curExp + tryWeights <- curWeights + tryDist <- oldInter + try.ess <- ess + print("OIOIOI") + } else{ + doResample <- FALSE + } + + out <- list(newExp=tryExp, weights=tryWeights, interDist=tryDist, ess=try.ess, doResample = doResample) return(out) } @@ -569,37 +516,69 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI # Auxiliary function for particle mutation -mutate <- function(setup, particles, proposalGenerator, posteriorValues, importanceDensity, method, steps, proposalScale, adaptive = TRUE){ +mutate <- function(setup, particles, proposalGenerator, posteriorValues, importanceDensity, method, steps, proposalScale, adaptive = TRUE, b){ + print(c("method in mutate", method)) + print(c("b",b)) if(is.vector(particles)){particles = matrix(particles, ncol = 1)} acceptance <- vector("numeric", length=steps) importanceValues <- importanceDensity(particles) - if(method=="Metropolis"){ + #if(method=="Metropolis"){ if(adaptive){ proposalGenerator = updateProposalGenerator(proposalGenerator, particles) } for(j in 1:steps){ - particlesProposals = proposalGenerator$returnProposalMatrix(particles, scale = proposalScale) + + if(method=="Metropolis"){ + particlesProposals = proposalGenerator$returnProposalMatrix(particles, scale = proposalScale) + + + } else if(method=="DE"){ + #particleDiffs <- vector("numeric", particleSize) + #randomTerms <- as.data.frame(matrix(ncol=ncol(particles), nrow = nrow(particles))) + particlesProposals <- particles + + for(part in 1:nrow(particles)){ + particlesOld <- particles + particleDiff <- rep(0,ncol(particles)) + newParts <- rep(0L,2) + + # Sample 2 other particles + while(all(particleDiff==0) | part %in% newParts){ + # 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 <- particlesOld[newParts[1],] - particlesOld[newParts[2],] + } + numPar <- ncol(particles) + randVector <- runif(numPar,-b,b) + #print(c("part", part)) + #print(c("Particle 1", head(particlesOld[newParts[1],]))) + #print(c("Particle 2", head(particlesOld[newParts[2],]))) + #print(c("particleDiff", head(particleDiff))) + particlesProposals[part,] <- particles[part,] + (particleDiff * proposalScale) + randVector + } + } + + print(head(particlesProposals)) + proposalPosteriors <- setup$posterior$density(particlesProposals) proposalImportance <- importanceDensity(particlesProposals) - #jumpProb <- exp(setup$posterior$density(particlesProposals) - likelihoodValues[sel])^(i/iterations) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) - #jumpProb <- exp(setup$posterior$density(particlesProposals) - setup$posterior$density(particles)) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) + print(c("proposalPosteriors", head(proposalPosteriors))) + print(c("posteriorValues", head(posteriorValues))) + jumpProb <- exp(proposalPosteriors - posteriorValues) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) - #jumpProb <- exp(setup$posterior$density(particlesProposals) - posteriorValues) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) - - - print(c("particlesProposals", head(particlesProposals))) - print(c("particles", head(particles))) + accepted <- jumpProb > runif(length(jumpProb), 0 ,1) - #rejectionRate = rejectionRate + sum(accepted) particles[accepted, ] = particlesProposals[accepted, ] posteriorValues[accepted] <- proposalPosteriors[accepted] importanceValues[accepted] <- proposalImportance[accepted] acceptance[j] <- sum(accepted)/nrow(particles) - #info$resamplingAcceptance[(icount),j] <- sum(accepted)/particleSize } - } + #} out <- list(particles=particles, posteriorValues = posteriorValues, importanceValues = importanceValues, acceptance = acceptance) } From 804e1954d87da750774fc041ecd91293b7c73813 Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Tue, 28 Aug 2018 10:28:12 +0200 Subject: [PATCH 14/42] Some unit tests for SMC --- BayesianTools/tests/testthat/test-SMC.R | 63 +++++++++++++++++++++++-- 1 file changed, 60 insertions(+), 3 deletions(-) diff --git a/BayesianTools/tests/testthat/test-SMC.R b/BayesianTools/tests/testthat/test-SMC.R index 856e1b7..54490bf 100644 --- a/BayesianTools/tests/testthat/test-SMC.R +++ b/BayesianTools/tests/testthat/test-SMC.R @@ -1,8 +1,15 @@ -skip_on_cran() +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) -context("Test SMC and utility functions") +library(testthat) +library(BayesianTools) +test <- "exact" +#skip_on_cran() + + + +set.seed(3) @@ -23,4 +30,54 @@ test_that("resampling works", { 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) + importanceDensity <- function(x){return(dunif(x, min=-10, max=10, log = TRUE))} + + posteriorValues <- setup$posterior$density(particles) + + # Initialize proposal Generator + proposalGenerator <- createProposalGenerator() + +}) \ No newline at end of file From 648e9496548dc53e3fd0a8d3bcd7029b474095c5 Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Tue, 28 Aug 2018 14:20:00 +0200 Subject: [PATCH 15/42] Add more diagnostic outputs to SMC --- BayesianTools/R/SMC.R | 129 +++++++----------------------------------- 1 file changed, 20 insertions(+), 109 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index 34e3477..e09ae8f 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -43,7 +43,8 @@ smcSampler <- function(bayesianSetup, pars.upper=NULL, mutate.method ="Metropolis", # TODO document b=1e-04, # TODO document - diagnostics = NULL){ + diagnostics = NULL, + reference=NULL){ ########### SETUP STEPS ######################## @@ -54,7 +55,8 @@ smcSampler <- function(bayesianSetup, # are made large enough that they will probably not need to be grown (growing arrays can be slow). 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$ess.vec <- info$exponents <- info$diagnostics <- vector("numeric", 10000) + diag.end <- vector("numeric", lastMutateSteps) setup <- checkBayesianSetup(bayesianSetup) @@ -159,75 +161,6 @@ smcSampler <- function(bayesianSetup, # 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 - #if (resampling == T){ - - - ####------------- Replaced by mutate function - Start ------------------- - # if(mutate=="Metropolis"){ - # - # if (adaptive == T){ - # proposalGenerator = updateProposalGenerator(proposalGenerator, particles) - # } - # - # for (j in 1:resamplingSteps){ - # particlesProposals = proposalGenerator$returnProposalMatrix(particles, scale = proposalScale) - # proposalPosteriors <- setup$posterior$density(particlesProposals) - # proposalImportance <- importanceDensity(particlesProposals) - # - # #jumpProb <- exp(setup$posterior$density(particlesProposals) - likelihoodValues[sel])^(i/iterations) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) - # #jumpProb <- exp(setup$posterior$density(particlesProposals) - setup$posterior$density(particles)) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) - # jumpProb <- exp(proposalPosteriors - posteriorValues) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) - # #jumpProb <- exp(setup$posterior$density(particlesProposals) - posteriorValues) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) - # - # - # print(c("particlesProposals", head(particlesProposals))) - # print(c("particles", head(particles))) - # accepted <- jumpProb > runif(length(jumpProb), 0 ,1) - # rejectionRate = rejectionRate + sum(accepted) - # particles[accepted, ] = particlesProposals[accepted, ] - # posteriorValues[accepted] <- proposalPosteriors[accepted] - # importanceValues[accepted] <- proposalImportance[accepted] - # info$resamplingAcceptance[(icount),j] <- sum(accepted)/particleSize - # } - # } else if(mutate=="DE"){ - # # Differential evolution - # # For now (test): basic DE algorithm - # - # for(j in 1:resamplingSteps){ - # - # # Loop over particles - # for(part in 1:particleSize){ - # particlesOld <- particles - # particleDiff <- rep(0,numPar) - # newParts <- rep(0L,2) - # - # # Sample 2 other particles - # while(all(particleDiff==0) | part %in% newParts){ - # # 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=particleSize,size=2) - # particleDiff <- particlesOld[newParts[1],] - particlesOld[newParts[2],] - # } - # - # particleProp <- particlesOld[part,] + (particleDiff * proposalScale) + runif(numPar,-b,b) - # jumpProb <- exp(setup$posterior$density(particleProp) - likelihoodValues[part]) * exp(setup$prior$density(particleProp) - setup$prior$density(particles[part,])) - # - # if(jumpProb > runif(length(jumpProb), 0 ,1)){ - # particles[part,] <- particleProp - # if(is.na(info$resamplingAcceptance[(icount-1),j])){ - # info$resamplingAcceptance[(icount-1),j] <- 1/particleSize - # } else{ - # info$resamplingAcceptance[(icount-1),j] <- info$resamplingAcceptance[(icount-1),j] + 1/particleSize - # } - # } else{ - # particles[part,] <- particlesOld[part,] - # } - # } - # } - # } - ####------------- Replaced by mutate function - End ------------------- 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 @@ -269,28 +202,6 @@ smcSampler <- function(bayesianSetup, print(c("interDist", head(interDist))) print(c("oldInter", head(oldInter))) - # Calculate new weights - #weights <- weights + (interDist - oldInter) - # Normalize (log-)weights so that the sum (of non-logs) equals 1 - #weights <- weights - BayesianTools:::logSumExp(weights) - - - # It may occur that the difference between the new and old weights of a particle is so great that an overflow occurs. - # In this case, the algorithm reduces the difference between the current and next intermediary distribution (i.e. the difference - # between the current and next value of the exponent). As a result, one more iteration is added. - # This is the same that happens when the effective sample size becomes *very* low. - #if(is.infinite(sum(sort(exp(weights))))){ - - # if(any(is.infinite(weights))){ - # print("Infinite weights sum; retrying with new intermediary distribution") - # expdiff <- (exponents[icount] - exponents[(icount-1)]) * 0.5 - # newExp <- exponents[(icount-1)] + expdiff - # exponents <- append(exponents,newExp,after=(icount-1)) - # weights <- oldweights - # info$res.of[icount] <- TRUE - # next - # } - print("-------") print(c("icount", icount)) @@ -345,25 +256,33 @@ smcSampler <- function(bayesianSetup, icount <- icount + 1 oldExp <- curExp - if(!is.null(diagnostics)) info$diagnostics[[icount-1]] <- diagnostics(particles) + if(!is.null(diagnostics)) info$diagnostics[[icount-1]] <- diagnostics(particles, reference) } # 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)) - # Last mutation, to increase diversity between particles - mutate.out <- mutate(setup = setup, particles = particles, proposalGenerator = proposalGenerator, posteriorValues = posteriorValues, importanceDensity = importanceDensity, method = mutate.method, steps = lastMutateSteps, proposalScale = proposalScale, adaptive = adaptive, b=b) - particles <- mutate.out$particles - posteriorValues <- mutate.out$posteriorValues - importanceValues <- mutate.out$importanceValues - lastAccept <- mutate.out$acceptance + 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) + particles <- mutate.out$particles + posteriorValues <- mutate.out$posteriorValues + importanceValues <- mutate.out$importanceValues + lastAccept <- mutate.out$acceptance + if(!is.null(diagnostics)) diag.end[mutateStep] <- diagnostics(particles, reference) + } info$rejectionRate = rejectionRate / (iterations * resamplingSteps) + # 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. 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$ess.vec <- info$ess.vec[1:(icount-1)] info$lastAccept <- lastAccept @@ -469,13 +388,7 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI while(abs(try.ess-target.ess) > tol & tryExp < 1 & abs(b-a) > 1e-10){ tryExp <- (a+b) * 0.5 - - # print(c("ess", ess)) - # print(c("target.ess", target.ess)) - # print(c("try.ess", try.ess)) - # print(c("curExp", curExp)) - # print(c("tryExp", tryExp)) - # print(c("a,b", a, b)) + tryDist <- tryExp * posteriorValues + (1-tryExp) * importanceValues @@ -485,8 +398,6 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI try.ess <- 1 / sum(exp(2 * tryWeights)) - #print(c("new try.ess", try.ess)) - 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 From 5960ac5fe8b6aa73e2c27bc76551d971240df2af Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Tue, 28 Aug 2018 16:01:24 +0200 Subject: [PATCH 16/42] Combine resample and mutate steps Mutation is only done after resampling --- BayesianTools/R/SMC.R | 46 ++++++++++--------------------------------- 1 file changed, 10 insertions(+), 36 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index e09ae8f..c52691c 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -1,5 +1,5 @@ #' SMC sampler -#' @author Florian Hartig +#' @author Florian Hartig, Matthias Speich #' @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. @@ -162,13 +162,11 @@ smcSampler <- function(bayesianSetup, # 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 - - + # 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 @@ -197,15 +195,6 @@ smcSampler <- function(bayesianSetup, curExp <- 1 } - print(c("weights U", head(weights))) - print(c("weights sum", BayesianTools:::logSumExp(weights))) - print(c("interDist", head(interDist))) - print(c("oldInter", head(oldInter))) - - print("-------") - print(c("icount", icount)) - - ######################################################### # Re?sampling step @@ -221,10 +210,6 @@ smcSampler <- function(bayesianSetup, # 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) - - print(c("ess", ess)) - print(c("ess.limit", ess.limit)) - print(c("curExp", curExp)) if(ess < ess.limit | icount == (length(exponents) - lastResample) | doResample){ @@ -414,7 +399,6 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI tryWeights <- curWeights tryDist <- oldInter try.ess <- ess - print("OIOIOI") } else{ doResample <- FALSE } @@ -427,14 +411,13 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI # Auxiliary function for particle mutation -mutate <- function(setup, particles, proposalGenerator, posteriorValues, importanceDensity, method, steps, proposalScale, adaptive = TRUE, b){ +mutate <- function(setup, particles, proposalGenerator, posteriorValues, importanceDensity, method, steps, proposalScale, adaptive = TRUE, b=1E-04){ print(c("method in mutate", method)) print(c("b",b)) if(is.vector(particles)){particles = matrix(particles, ncol = 1)} acceptance <- vector("numeric", length=steps) importanceValues <- importanceDensity(particles) - #if(method=="Metropolis"){ if(adaptive){ proposalGenerator = updateProposalGenerator(proposalGenerator, particles) } @@ -446,8 +429,7 @@ mutate <- function(setup, particles, proposalGenerator, posteriorValues, importa } else if(method=="DE"){ - #particleDiffs <- vector("numeric", particleSize) - #randomTerms <- as.data.frame(matrix(ncol=ncol(particles), nrow = nrow(particles))) + particlesProposals <- particles for(part in 1:nrow(particles)){ @@ -466,22 +448,14 @@ mutate <- function(setup, particles, proposalGenerator, posteriorValues, importa } numPar <- ncol(particles) randVector <- runif(numPar,-b,b) - #print(c("part", part)) - #print(c("Particle 1", head(particlesOld[newParts[1],]))) - #print(c("Particle 2", head(particlesOld[newParts[2],]))) - #print(c("particleDiff", head(particleDiff))) particlesProposals[part,] <- particles[part,] + (particleDiff * proposalScale) + randVector } } - print(head(particlesProposals)) - + proposalPosteriors <- setup$posterior$density(particlesProposals) proposalImportance <- importanceDensity(particlesProposals) - - print(c("proposalPosteriors", head(proposalPosteriors))) - print(c("posteriorValues", head(posteriorValues))) - + jumpProb <- exp(proposalPosteriors - posteriorValues) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) accepted <- jumpProb > runif(length(jumpProb), 0 ,1) From 175a8937edf3824ccd89f4479ae8a3a3f0afeab8 Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Wed, 29 Aug 2018 11:16:15 +0200 Subject: [PATCH 17/42] Allow SIR without mutation --- BayesianTools/R/SMC.R | 81 +++++++++++++++++++++++-------------------- 1 file changed, 43 insertions(+), 38 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index c52691c..2175cde 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -48,12 +48,12 @@ smcSampler <- function(bayesianSetup, ########### SETUP STEPS ######################## - if(resamplingSteps < 1) stop("SMC error, resamplingSteps can't be < 1") + #if(resamplingSteps < 1) stop("SMC error, resamplingSteps can't be < 1") info = list() # 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). - 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 + 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 <- info$diagnostics <- vector("numeric", 10000) diag.end <- vector("numeric", lastMutateSteps) @@ -73,7 +73,13 @@ smcSampler <- function(bayesianSetup, } if(class(initialParticles) == "list"){ importanceDensity = initialParticles$density - initialParticles = as.matrix(initialParticles$particles,ncol=1) + 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") @@ -145,7 +151,7 @@ smcSampler <- function(bayesianSetup, #while(curExp <= 1 & !lastIteration){ # For fully adaptive algorithm # if(curExp == 1){lastIteration <- TRUE} while(curExp < 1){ - + if (numPar == 1) particles = matrix(particles, ncol = 1) @@ -183,12 +189,15 @@ smcSampler <- function(bayesianSetup, # Create intermediary distribution ess <- 1 / sum(exp(2 * weights)) oldInter <- oldExp * posteriorValues + (1-oldExp) * importanceValues + inter.out <- beta.search(ess=ess, target.ess = (ess * ess.factor), 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){ @@ -197,16 +206,7 @@ smcSampler <- function(bayesianSetup, ######################################################### # Re?sampling step - - # Option 1) Always resample (default BT) - # Option 2) Keep weights and resample if ESS < something - # => for now, option 1 is activated if the argument ess.limit is set to a value >= the number of particles - - # Question: last step complete resampling, so that all particles are the same? - - if(is.infinite(ess)){print(weights)} - # 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) @@ -216,8 +216,9 @@ smcSampler <- function(bayesianSetup, oldExp <- curExp sel = resample(weights, method = sampling) - + particles = particles[sel,] + posteriorValues <- posteriorValues[sel] importanceValues <- importanceValues[sel] info$survivingParticles[icount] = length(unique(sel)) @@ -227,12 +228,15 @@ smcSampler <- function(bayesianSetup, # Normalize (log-)weights so that the sum (of non-logs) equals 1 weights <- weights - BayesianTools:::logSumExp(weights) - print("-- Resampling") - 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 + + 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) + particles <- mutate.out$particles + posteriorValues <- mutate.out$posteriorValues + importanceValues <- mutate.out$importanceValues + info$resamplingAcceptance[(icount),] <- mutate.out$acceptance + } + } oldweights <- weights @@ -244,32 +248,35 @@ smcSampler <- function(bayesianSetup, 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)) - 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) - particles <- mutate.out$particles - posteriorValues <- mutate.out$posteriorValues - importanceValues <- mutate.out$importanceValues - lastAccept <- mutate.out$acceptance - if(!is.null(diagnostics)) diag.end[mutateStep] <- diagnostics(particles, reference) + 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) + particles <- mutate.out$particles + posteriorValues <- mutate.out$posteriorValues + importanceValues <- mutate.out$importanceValues + lastAccept <- mutate.out$acceptance + if(!is.null(diagnostics)) diag.end[mutateStep] <- diagnostics(particles, reference) + } } - info$rejectionRate = rejectionRate / (iterations * resamplingSteps) # 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. - info$resamplingAcceptance <- info$resamplingAcceptance[1:(icount-1),] + 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$ess.vec <- info$ess.vec[1:(icount-1)] - info$lastAccept <- lastAccept + if(resamplingSteps > 0) info$lastAccept <- lastAccept settings = list(initialParticles = initialParticles, proposalGenerator = proposalGenerator) @@ -366,10 +373,7 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI b <- 1 try.ess <- target.ess + (100*tol) # Dummy initial value that is sure to be different from target value tryExp <- 0 - - print("*-*-*_*-*-*") - print(c("ess", ess)) - print(c("curWeights", head(curWeights))) + while(abs(try.ess-target.ess) > tol & tryExp < 1 & abs(b-a) > 1e-10){ tryExp <- (a+b) * 0.5 @@ -379,8 +383,8 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI tryWeights <- curWeights + (tryDist - oldInter) # Normalize (log-)weights so that the sum (of non-logs) equals 1 + tryWeights <- tryWeights - BayesianTools:::logSumExp(tryWeights) - try.ess <- 1 / sum(exp(2 * tryWeights)) if(try.ess - target.ess > tol & !any(is.infinite(tryWeights))){ @@ -392,6 +396,8 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI b <- tryExp } } + + if(any(is.infinite(tryWeights)) | is.infinite(try.ess)){ doResample <- TRUE @@ -412,8 +418,7 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI mutate <- function(setup, particles, proposalGenerator, posteriorValues, importanceDensity, method, steps, proposalScale, adaptive = TRUE, b=1E-04){ - print(c("method in mutate", method)) - print(c("b",b)) + if(is.vector(particles)){particles = matrix(particles, ncol = 1)} acceptance <- vector("numeric", length=steps) importanceValues <- importanceDensity(particles) From f5710a1ca5cc7e493c89fd2d5b152151dadfe1c3 Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Wed, 29 Aug 2018 11:16:37 +0200 Subject: [PATCH 18/42] Add more unit tests for SMC --- BayesianTools/tests/testthat/test-SMC.R | 72 +++++++++++++++++++++++-- 1 file changed, 68 insertions(+), 4 deletions(-) diff --git a/BayesianTools/tests/testthat/test-SMC.R b/BayesianTools/tests/testthat/test-SMC.R index 54490bf..77ff58c 100644 --- a/BayesianTools/tests/testthat/test-SMC.R +++ b/BayesianTools/tests/testthat/test-SMC.R @@ -2,6 +2,7 @@ context("Test SMC and utility functions") # Must be at the top of the file (http library(testthat) library(BayesianTools) +library(Matching) test <- "exact" @@ -9,7 +10,7 @@ test <- "exact" -set.seed(3) +set.seed(123) @@ -73,11 +74,74 @@ test_that("mutation works", { 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))} - posteriorValues <- setup$posterior$density(particles) + ### Testing *Metropolis* implementation + posteriorValues <- setup$posterior$density(particles) # Initialize proposal Generator - proposalGenerator <- createProposalGenerator() + 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) -}) \ No newline at end of file + # 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) +}) From 3465ae18f3273c18bcf16a26b297f074b24a9b75 Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Mon, 8 Oct 2018 16:46:14 +0200 Subject: [PATCH 19/42] Lower tolerance criterion to avoid crash of Bhattacharyya distance function Tolerance parameter (passed to solve) is now set to the smallest number that can be represented on the machine. --- BayesianTools/R/Distance.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/BayesianTools/R/Distance.R b/BayesianTools/R/Distance.R index c644e82..b8f7897 100644 --- a/BayesianTools/R/Distance.R +++ b/BayesianTools/R/Distance.R @@ -30,7 +30,8 @@ getSampleDistance <- function(sample1, sample2, type = "KL"){ # The following code is copied from package fpc aggregatesigma <- (Sigma1+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=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 out <- d1+d2 From 93da4169b988659695caf9e45ae994abfc8e38ae Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Mon, 8 Oct 2018 16:49:23 +0200 Subject: [PATCH 20/42] Improve proposal function of MCMC kernel - Automatic selection of scaling factors - Snooker update --- BayesianTools/R/SMC.R | 326 ++++++++++++++++++++++++++++++++++++++---- 1 file changed, 300 insertions(+), 26 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index 2175cde..fdd2a0f 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -48,6 +48,10 @@ smcSampler <- function(bayesianSetup, ########### SETUP STEPS ######################## + # 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() @@ -56,7 +60,7 @@ smcSampler <- function(bayesianSetup, 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 <- info$diagnostics <- vector("numeric", 10000) - diag.end <- vector("numeric", lastMutateSteps) + diag.end <- lastAccept <- vector("numeric", lastMutateSteps) setup <- checkBayesianSetup(bayesianSetup) @@ -99,6 +103,10 @@ smcSampler <- function(bayesianSetup, } usedUp = 0 + + #### TEST + proposalScale <- runif(nrow(particles), min = 1e-06, max = 1) + Z <- particles ########### WEIGHT SETUP ######################## @@ -128,7 +136,9 @@ smcSampler <- function(bayesianSetup, # Initial importance distribution importanceValues <- importanceDensity(particles) + #print("Calculating initial posteriors") posteriorValues <- setup$posterior$density(particles) + #print("Initial posteriors done") #oldInter <- importanceValues oldExp <- 0 @@ -141,6 +151,10 @@ smcSampler <- function(bayesianSetup, 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) + ########################################################## # Loop starts here # @@ -151,10 +165,15 @@ smcSampler <- function(bayesianSetup, #while(curExp <= 1 & !lastIteration){ # For fully adaptive algorithm # if(curExp == 1){lastIteration <- TRUE} while(curExp < 1){ + + # 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) - # Using a while loop instead of for because in the adaptive algorithm, the number of iterations may increase ########################################################### @@ -196,6 +215,7 @@ smcSampler <- function(bayesianSetup, interDist <- inter.out$interDist ess <- inter.out$ess doResample <- inter.out$doResample + info$ess.vec[icount] <- ess @@ -211,13 +231,16 @@ smcSampler <- function(bayesianSetup, # 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 | icount == (length(exponents) - lastResample) | doResample){ + if(ess < ess.limit | 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] @@ -229,16 +252,60 @@ smcSampler <- function(bayesianSetup, 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) + 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 = 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", 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 @@ -255,19 +322,31 @@ smcSampler <- function(bayesianSetup, 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) + 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 <- mutate.out$acceptance + 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. @@ -277,6 +356,7 @@ smcSampler <- function(bayesianSetup, info$diagnostics <- c(info$diagnostics[1:(icount-1)], 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) @@ -293,6 +373,10 @@ smcSampler <- function(bayesianSetup, } +###################################################### +# END of SMC function +###################################################### + ###################################################### # Auxiliary functions for resampling @@ -333,7 +417,7 @@ systematicResampling <- function(weights){ U[i] <- ((i-1)/n.parts) + u - new.parts.sort[i] <- tail(which(cumu.weights <= U[i]),1) + 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]) } @@ -417,12 +501,27 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI # Auxiliary function for particle mutation -mutate <- function(setup, particles, proposalGenerator, posteriorValues, importanceDensity, method, steps, proposalScale, adaptive = TRUE, b=1E-04){ +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) } @@ -431,44 +530,219 @@ mutate <- function(setup, particles, proposalGenerator, posteriorValues, importa if(method=="Metropolis"){ particlesProposals = proposalGenerator$returnProposalMatrix(particles, scale = proposalScale) - - } else if(method=="DE"){ particlesProposals <- particles for(part in 1:nrow(particles)){ - particlesOld <- particles - particleDiff <- rep(0,ncol(particles)) - newParts <- rep(0L,2) - # Sample 2 other particles - while(all(particleDiff==0) | part %in% newParts){ - # 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 <- particlesOld[newParts[1],] - particlesOld[newParts[2],] + 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])) + #} } - numPar <- ncol(particles) - randVector <- runif(numPar,-b,b) - particlesProposals[part,] <- particles[part,] + (particleDiff * proposalScale) + randVector } } 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) * exp(setup$prior$density(particlesProposals) - setup$prior$density(particles)) + jumpProb <- exp(proposalPosteriors - posteriorValues + r_extra) accepted <- jumpProb > runif(length(jumpProb), 0 ,1) 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)) + } #} - out <- list(particles=particles, posteriorValues = posteriorValues, importanceValues = importanceValues, acceptance = acceptance) + 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 <- function(particle.row, min.pars, max.pars){ + scaled.particles <- (particle.row - min.pars) / (max.pars - min.pars) + return(scaled.particles) } + +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 +# +# } \ No newline at end of file From a3535f5f410badc875d811a187d2c20b7ea2e827 Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Tue, 16 Oct 2018 16:57:04 +0200 Subject: [PATCH 21/42] Add possibility to get two output terms for Bhattacharyya distance The two output terms represent the contributions of means and covariance to total distance --- BayesianTools/R/Distance.R | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/BayesianTools/R/Distance.R b/BayesianTools/R/Distance.R index b8f7897..4bba935 100644 --- a/BayesianTools/R/Distance.R +++ b/BayesianTools/R/Distance.R @@ -21,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) @@ -34,7 +34,11 @@ getSampleDistance <- function(sample1, sample2, type = "KL"){ 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 From 25981c43d80f6ec7385cfb38a3c91403a38158db Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Tue, 16 Oct 2018 16:59:18 +0200 Subject: [PATCH 22/42] Allow 2 output terms for diagnostic function --- BayesianTools/R/SMC.R | 46 +++++++++++++++++++++++-------------------- 1 file changed, 25 insertions(+), 21 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index fdd2a0f..883bbab 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -59,8 +59,9 @@ smcSampler <- function(bayesianSetup, # 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 <- info$diagnostics <- vector("numeric", 10000) - diag.end <- lastAccept <- vector("numeric", lastMutateSteps) + info$ess.vec <- info$exponents <- vector("numeric", 10000) + info$diagnostics <- diag.end <- list() + lastAccept <- vector("numeric", lastMutateSteps) setup <- checkBayesianSetup(bayesianSetup) @@ -105,7 +106,7 @@ smcSampler <- function(bayesianSetup, usedUp = 0 #### TEST - proposalScale <- runif(nrow(particles), min = 1e-06, max = 1) + #proposalScale <- runif(nrow(particles), min = 1e-06, max = 1) Z <- particles @@ -209,7 +210,8 @@ smcSampler <- function(bayesianSetup, ess <- 1 / sum(exp(2 * weights)) oldInter <- oldExp * posteriorValues + (1-oldExp) * importanceValues - 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 * ess.factor), posteriorValues = posteriorValues, importanceValues = importanceValues, oldInter = oldInter, curWeights = weights, curExp = curExp) + inter.out <- beta.search(ess=ess, target.ess = (ess.limit-1), posteriorValues = posteriorValues, importanceValues = importanceValues, oldInter = oldInter, curWeights = weights, curExp = curExp) curExp <- inter.out$newExp weights <- inter.out$weights interDist <- inter.out$interDist @@ -274,22 +276,23 @@ smcSampler <- function(bayesianSetup, proposalScale <- mutate.out$proposalScale } - plot(density(particles[,1]), xlim=c(0.314, 0.681), main = nrow(unique(particles))) + #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", nrow(unique(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)))) + # } } @@ -333,9 +336,9 @@ smcSampler <- function(bayesianSetup, posteriorValues <- mutate.out$posteriorValues importanceValues <- mutate.out$importanceValues lastAccept[mutateStep] <- mutate.out$acceptance - if(!is.null(diagnostics)) diag.end[mutateStep] <- diagnostics(particles, reference) + if(!is.null(diagnostics)) diag.end[[mutateStep]] <- diagnostics(particles, reference) - plot(density(particles[,1]), xlim=c(0.314, 0.681)) + #plot(density(particles[,1]), xlim=c(0.314, 0.681)) #Z <- rbind(Z, particles) } @@ -353,7 +356,8 @@ smcSampler <- function(bayesianSetup, 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[[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 @@ -691,7 +695,7 @@ mutate <- function(setup, particles, proposalGenerator, posteriorValues, importa } - #print(head(scale.factors, 100)) + print(head(scale.factors, 100)) #plot(density(scale.factors), xlim=c(0,1)) } From c90cff2659224031c9dacd42e46dbffcf624efee Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Tue, 16 Oct 2018 18:35:59 +0200 Subject: [PATCH 23/42] Suppress print --- BayesianTools/R/SMC.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index 883bbab..6645838 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -695,7 +695,7 @@ mutate <- function(setup, particles, proposalGenerator, posteriorValues, importa } - print(head(scale.factors, 100)) + #print(head(scale.factors, 100)) #plot(density(scale.factors), xlim=c(0,1)) } From 72ec8efcb012a8a75853aca32f772118a81603f2 Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Mon, 5 Nov 2018 16:14:14 +0100 Subject: [PATCH 24/42] Fix code for Distance (type "D") --- BayesianTools/R/Distance.R | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/BayesianTools/R/Distance.R b/BayesianTools/R/Distance.R index 4bba935..e36b747 100644 --- a/BayesianTools/R/Distance.R +++ b/BayesianTools/R/Distance.R @@ -32,8 +32,7 @@ getSampleDistance <- function(sample1, sample2, type = "KL"){ aggregatesigma <- (Sigma1+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 + 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"){ @@ -50,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 ) From 0245c0a2cb39d198f743d1464e69a20d8fa2ad31 Mon Sep 17 00:00:00 2001 From: Florian Hartig Date: Tue, 13 Nov 2018 13:30:15 +0100 Subject: [PATCH 25/42] changed runtimeto info slot --- BayesianTools/R/classMcmcSampler.R | 2 +- BayesianTools/R/classMcmcSamplerList.R | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/BayesianTools/R/classMcmcSampler.R b/BayesianTools/R/classMcmcSampler.R index 8395fb7..59bca43 100644 --- a/BayesianTools/R/classMcmcSampler.R +++ b/BayesianTools/R/classMcmcSampler.R @@ -131,7 +131,7 @@ summary.mcmcSampler <- function(object, ...){ 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, ...) diff --git a/BayesianTools/R/classMcmcSamplerList.R b/BayesianTools/R/classMcmcSamplerList.R index f0a8dbe..859ccd9 100644 --- a/BayesianTools/R/classMcmcSamplerList.R +++ b/BayesianTools/R/classMcmcSamplerList.R @@ -33,7 +33,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) From 57eabae695aa4c7e6cc0b9ca398a0d340c34e6be Mon Sep 17 00:00:00 2001 From: Florian Hartig Date: Tue, 13 Nov 2018 13:30:43 +0100 Subject: [PATCH 26/42] update to sample distance help --- BayesianTools/inst/examples/getSampleDistanceHelp.R | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/BayesianTools/inst/examples/getSampleDistanceHelp.R b/BayesianTools/inst/examples/getSampleDistanceHelp.R index b78c721..647200a 100644 --- a/BayesianTools/inst/examples/getSampleDistanceHelp.R +++ b/BayesianTools/inst/examples/getSampleDistanceHelp.R @@ -12,3 +12,11 @@ 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") From 9a3f2917d0c2b45ae4e99e77e9b8373f80263139 Mon Sep 17 00:00:00 2001 From: Florian Hartig Date: Tue, 13 Nov 2018 13:33:06 +0100 Subject: [PATCH 27/42] write runtime info to info slot --- BayesianTools/R/mcmcRun.R | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/BayesianTools/R/mcmcRun.R b/BayesianTools/R/mcmcRun.R index de717ed..3ca9671 100644 --- a/BayesianTools/R/mcmcRun.R +++ b/BayesianTools/R/mcmcRun.R @@ -236,9 +236,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) } @@ -448,9 +450,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)){ From 227d047937487571fb801073c77baf3a22665437 Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Tue, 13 Nov 2018 14:50:55 +0100 Subject: [PATCH 28/42] Temporary changes for debugging --- .gitignore | 3 ++- BayesianTools/R/SMC.R | 19 +++++++++++++++++++ 2 files changed, 21 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index ff72220..fbb14fa 100644 --- a/.gitignore +++ b/.gitignore @@ -49,4 +49,5 @@ examples/16-TG14/TG14-v2.2_files/figure-html BayesianTools/vignettes/BayesianTools_cache BayesianTools/src/vsem.o-6bf934a9 -Questions/ \ No newline at end of file +Questions/ +BayesianTools/R/test_basfor.r diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index 6645838..b18fbfa 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -471,6 +471,25 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI 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)) From 7f5c81b8f4583186b72779a7cdd3025410b080f6 Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Tue, 13 Nov 2018 17:57:00 +0100 Subject: [PATCH 29/42] Temporary outputs for debugging --- BayesianTools/R/SMC.R | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index b18fbfa..aaae1a3 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -649,6 +649,18 @@ mutate <- function(setup, particles, proposalGenerator, posteriorValues, importa 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] From a5ec62fe91eb5b36f4d6ca75ca5288f84464044a Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Wed, 14 Nov 2018 15:06:04 +0100 Subject: [PATCH 30/42] Temporary debug output --- BayesianTools/R/SMC.R | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index aaae1a3..7348888 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -167,6 +167,14 @@ smcSampler <- function(bayesianSetup, # if(curExp == 1){lastIteration <- TRUE} while(curExp < 1){ + ### TEMPORARY DEBUG + if(sum(is.nan(importanceValues)) > 1){ + print("NaN in importanceValues") + print(c("indices", which(is.nan(importanceValues)))) + print("Particles with NaN:") + print(particles[is.nan(importanceValues),]) + } + # 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) From fcd31972c41a57a5fa19e8281b83b9eab9226da0 Mon Sep 17 00:00:00 2001 From: Florian Hartig Date: Wed, 14 Nov 2018 16:32:05 +0100 Subject: [PATCH 31/42] distance help update --- .../inst/examples/getSampleDistanceHelp.R | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/BayesianTools/inst/examples/getSampleDistanceHelp.R b/BayesianTools/inst/examples/getSampleDistanceHelp.R index 647200a..6f06340 100644 --- a/BayesianTools/inst/examples/getSampleDistanceHelp.R +++ b/BayesianTools/inst/examples/getSampleDistanceHelp.R @@ -20,3 +20,18 @@ 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") From 863464652a16a0fd516507f8703af34b50222e0a Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Tue, 20 Nov 2018 16:06:42 +0100 Subject: [PATCH 32/42] Temporary debug output --- BayesianTools/R/SMC.R | 27 +++++++++++++++++++++------ 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index 7348888..36d98a2 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -168,11 +168,25 @@ smcSampler <- function(bayesianSetup, while(curExp < 1){ ### TEMPORARY DEBUG - if(sum(is.nan(importanceValues)) > 1){ - print("NaN in importanceValues") - print(c("indices", which(is.nan(importanceValues)))) - print("Particles with NaN:") - print(particles[is.nan(importanceValues),]) + 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)) > 1){ + print(c("NaN indices", which(is.nan(posteriorValues)))) + print("Particles with NaN:") + print(particles[is.nan(posteriorValues),]) + } + + if(sum(is.na(posteriorValues)) > 1){ + print(c("NA indices", which(is.na(posteriorValues)))) + print("Particles with NA:") + print(particles[is.na(posteriorValues),]) + } + + if(sum(is.infinite(posteriorValues)) > 1){ + 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) @@ -461,7 +475,7 @@ resample <- function(weights, method = "multinomial"){ 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. (2001), Scand J Statist, doi: 10.1111/j.1467-9469.2010.00723.x + # 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 @@ -497,6 +511,7 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI print(posteriorValues) print("importanceValues") print(importanceValues) + } tryWeights <- tryWeights - BayesianTools:::logSumExp(tryWeights) From 5e4165f99153edc7c63e55d509bc6c781ea1a9ef Mon Sep 17 00:00:00 2001 From: Matthias Speich Date: Tue, 20 Nov 2018 16:29:45 +0100 Subject: [PATCH 33/42] Fix temporary debug output --- BayesianTools/R/SMC.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index 36d98a2..8f07779 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -170,19 +170,19 @@ smcSampler <- function(bayesianSetup, ### 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)) > 1){ + 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)) > 1){ + 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)) > 1){ + if(sum(is.infinite(posteriorValues)) > 0){ print(c("Infinite indices", which(is.infinite(posteriorValues)))) print("Particles with Infinite:") print(particles[is.infinite(posteriorValues),]) From 62377758c9627d55ae97f84700ccaf81af09e59b Mon Sep 17 00:00:00 2001 From: florianhartig Date: Mon, 19 Apr 2021 18:39:49 +0200 Subject: [PATCH 34/42] Updated help / documentation --- BayesianTools/R/SMC.R | 53 +++++++++++++++++++++++-------------------- 1 file changed, 29 insertions(+), 24 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index 89e2ff7..7d48edf 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -5,7 +5,8 @@ #' @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 if new particles should be created at each iteration -#' @param resamplingSteps how many resampling (MCMC) steps between the iterations +#' @param resamplingSteps how many resampling (MCMC) steps between the SMC iterations +#' @param lastMutateSteps how many resampling (MCMC) steps after the SMC iterations #' @param proposal optional proposal class #' @param exponents series of exponents to build the intermediate distributions #' @param adaptive should the covariance of the proposal be adapted during sampling @@ -14,9 +15,12 @@ #' @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 of effective sample size below which resampling is done. By default, the value is set to half the number of particles. To resample at each step, use a value >= the number of particles. +#' @param ess.factor #' @param lastResample Iteration (starting from the end) at which particle resampling is forced. To deactivate this, set to a value < 0 #' @param pars.lower Optional: vector contaiing the mimimum values of calibration parameters. Used to initialize the proposal function #' @param pars.upper Optional: vector contaiing the maximum values of calibration parameters. Used to initialize the proposal function +#' @param mutate.method MCMC method for the resampling. Either "Metropolis" or "DE" +#' @param b parameter for snooker update in DE #' @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. @@ -28,7 +32,7 @@ smcSampler <- function(bayesianSetup, iterations = 10, resampling = T, resamplingSteps = 2, - lastMutateSteps = 5, # TODO document + lastMutateSteps = 5, proposal = NULL, exponents = NULL, adaptive = T, @@ -46,7 +50,7 @@ smcSampler <- function(bayesianSetup, diagnostics = NULL, reference=NULL){ - ########### SETUP STEPS ######################## + # SETUP STEPS # Timing: normally done in mcmcRun.R. Timing is included in this function for the purpose of testing the SMC function. # Starting the clock @@ -110,7 +114,7 @@ smcSampler <- function(bayesianSetup, Z <- particles - ########### WEIGHT SETUP ######################## + # 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 @@ -157,10 +161,6 @@ smcSampler <- function(bayesianSetup, max.pars <- apply(particles,2,max) - ########################################################## - # Loop starts here # - ########################################################## - #for (i in 1:iterations){ # while(icount <= length(exponents)){ #while(curExp <= 1 & !lastIteration){ # For fully adaptive algorithm @@ -199,7 +199,6 @@ smcSampler <- function(bayesianSetup, # 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 @@ -248,7 +247,6 @@ smcSampler <- function(bayesianSetup, curExp <- 1 } - ######################################################### # Re?sampling step # Determine if resampling is necessary - if yes, resample with probability given by the weight @@ -399,12 +397,6 @@ smcSampler <- function(bayesianSetup, } -###################################################### -# END of SMC function -###################################################### - -###################################################### -# Auxiliary functions for resampling #' Residual Resampling #' @@ -452,7 +444,10 @@ systematicResampling <- function(weights){ } - +#' Auxiliary function for adaptive SMC algorithm +#' +#' +#' @keywords internal resample <- function(weights, method = "multinomial"){ particleSize = length(weights) @@ -469,10 +464,11 @@ resample <- function(weights, method = "multinomial"){ return(sel) } -###################################################### - -# Auxiliary function for adaptive algorithm +#' 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 @@ -543,10 +539,11 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI return(out) } -###################################################### -# Auxiliary function for particle mutation - +#' 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)} @@ -763,11 +760,19 @@ mutate <- function(setup, particles, proposalGenerator, posteriorValues, importa 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) @@ -797,7 +802,7 @@ rao.div <- function(scaled.particles){ } -###### + # Auxiliary function to calculate Metropolis ratio when log-likelihood values are very large # largeRatio <- function(a,b){ From 3fed011c8caf2e957325a9d9bcff71d6d400ebef Mon Sep 17 00:00:00 2001 From: florianhartig Date: Mon, 19 Apr 2021 19:28:23 +0200 Subject: [PATCH 35/42] help updates --- BayesianTools/DESCRIPTION | 2 +- BayesianTools/inst/examples/SMCHelp.R | 88 +--------------- BayesianTools/man/DREAM.Rd | 5 +- BayesianTools/man/beta.search.Rd | 21 ++++ BayesianTools/man/calibrationTest.Rd | 15 ++- BayesianTools/man/getSampleDistance.Rd | 23 ++++ BayesianTools/man/mutate.Rd | 26 +++++ BayesianTools/man/rao.div.Rd | 12 +++ BayesianTools/man/resample.Rd | 12 +++ BayesianTools/man/residualResampling.Rd | 12 +++ BayesianTools/man/scale.particles.Rd | 12 +++ BayesianTools/man/smcSampler.Rd | 134 +++++++++--------------- 12 files changed, 183 insertions(+), 179 deletions(-) create mode 100644 BayesianTools/man/beta.search.Rd create mode 100644 BayesianTools/man/mutate.Rd create mode 100644 BayesianTools/man/rao.div.Rd create mode 100644 BayesianTools/man/resample.Rd create mode 100644 BayesianTools/man/residualResampling.Rd create mode 100644 BayesianTools/man/scale.particles.Rd diff --git a/BayesianTools/DESCRIPTION b/BayesianTools/DESCRIPTION index e2e475f..32f7e6c 100644 --- a/BayesianTools/DESCRIPTION +++ b/BayesianTools/DESCRIPTION @@ -40,7 +40,7 @@ Suggests: testthat LinkingTo: Rcpp LazyData: true -RoxygenNote: 7.0.2.9000 +RoxygenNote: 7.1.1 URL: https://github.com/florianhartig/BayesianTools BugReports: https://github.com/florianhartig/BayesianTools/issues VignetteBuilder: knitr diff --git a/BayesianTools/inst/examples/SMCHelp.R b/BayesianTools/inst/examples/SMCHelp.R index 8270740..fcf4dc2 100644 --- a/BayesianTools/inst/examples/SMCHelp.R +++ b/BayesianTools/inst/examples/SMCHelp.R @@ -2,93 +2,13 @@ # 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), +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) - -} - - - +SMCout <- smcSampler(bayesianSetup) +marginalPlot(SMCout, type = "d", singlePanel = T) diff --git a/BayesianTools/man/DREAM.Rd b/BayesianTools/man/DREAM.Rd index 6aad1d6..7d7179a 100644 --- a/BayesianTools/man/DREAM.Rd +++ b/BayesianTools/man/DREAM.Rd @@ -7,8 +7,9 @@ DREAM( bayesianSetup, settings = list(iterations = 10000, nCR = 3, gamma = NULL, eps = 0, e = 0.05, - pCRupdate = TRUE, updateInterval = 10, burnin = 0, thin = 1, adaptation = 0.2, DEpairs - = 2, consoleUpdates = 10, startValue = NULL, currentChain = 1, message = TRUE) + pCRupdate = TRUE, updateInterval = 10, burnin = 0, thin = 1, adaptation = 0.2, + DEpairs = 2, consoleUpdates = 10, startValue = NULL, currentChain = 1, message = + TRUE) ) } \arguments{ 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/calibrationTest.Rd b/BayesianTools/man/calibrationTest.Rd index 75edcf3..e548e0f 100644 --- a/BayesianTools/man/calibrationTest.Rd +++ b/BayesianTools/man/calibrationTest.Rd @@ -7,14 +7,14 @@ calibrationTest(posteriorList, priorDraws, ...) } \arguments{ -\item{posteriorList}{a list with posterior samples. List items must be of a class that is supported by \code{\link{getSample}}} +\item{posteriorList}{a list with posterior samples. List items must be of a class that is supported by \code{\link{getSample}}. This includes BayesianTools objects, but also matrix and data.frame} \item{priorDraws}{a matrix with parameter values, drawn from the prior, that were used to simulate the data underlying the posteriorList. If colnames are provided, these will be used in the plots} \item{...}{arguments to be passed to \code{\link{getSample}}. Consider in particular the thinning option.} } \description{ -This function data averaged posterior +This function performs simulation-based calibration tests based on the idea that posteriors averaged over the prior should yield the prior. } \details{ The purpose of this function is to evaluate the results of a simulation-based calibration of an MCMC analysis. @@ -25,14 +25,19 @@ Briefly, the idea is to repeatedly 2. simulate new data based on these parameters, 3. calculate the posterior for these data -If the sampler and the likelihood are implemented correctly, the average of over all the posterior distribution should then again yield the prior (e.g. Cook et al., 2006). +If the sampler and the likelihood are implemented correctly, the average over all the posterior distribution should then again yield the prior (e.g. Cook et al., 2006). To test if this is the case, we implement the methods suggested by Talts et al., which is to calculate the rank statistics between the parameter draws and the posterior draws, which we then formally evaluate with a qq unif plot, and a ks.test -We speculate that a ks.test between the two distribution would likely give an identical result, but this is not noted in Talts et al. - +I speculate that a ks.test between the two distribution would likely give an identical result, but this is not noted in Talts et al. Cook, S. R., Gelman, A. and Rubin, D. B. (2006). Validation of Software for Bayesian Models Using Posterior Quantiles. J. Comput. Graph. Stat. 15 675-692. Talts, Sean, Michael Betancourt, Daniel Simpson, Aki Vehtari, and Andrew Gelman. "Validating Bayesian Inference Algorithms with Simulation-Based Calibration." arXiv preprint arXiv:1804.06788 (2018). } +\note{ +This function was implemented for the tests in Maliet, Odile, Florian Hartig, and Hélène Morlon. "A model with many small shifts for estimating species-specific diversification rates." Nature ecology & evolution 3.7 (2019): 1086-1092. The code linked with this paper provides a further example of its use. +} +\author{ + +} diff --git a/BayesianTools/man/getSampleDistance.Rd b/BayesianTools/man/getSampleDistance.Rd index 3d3d1fd..2dc346d 100644 --- a/BayesianTools/man/getSampleDistance.Rd +++ b/BayesianTools/man/getSampleDistance.Rd @@ -38,4 +38,27 @@ 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 31382f1..eb2b4f3 100644 --- a/BayesianTools/man/smcSampler.Rd +++ b/BayesianTools/man/smcSampler.Rd @@ -10,9 +10,23 @@ smcSampler( iterations = 10, resampling = T, 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 = NULL, + ess.factor = 0.95, + lastResample = 1, + pars.lower = NULL, + pars.upper = NULL, + mutate.method = "Metropolis", + b = 1e-04, + diagnostics = NULL, + reference = NULL ) } \arguments{ @@ -24,13 +38,39 @@ smcSampler( \item{resampling}{if new particles should be created at each iteration} -\item{resamplingSteps}{how many resampling (MCMC) steps between the iterations} +\item{resamplingSteps}{how many resampling (MCMC) steps between the SMC iterations} + +\item{lastMutateSteps}{how many resampling (MCMC) steps after the SMC iterations} \item{proposal}{optional proposal class} +\item{exponents}{series of exponents to build the intermediate distributions} + \item{adaptive}{should the covariance of the proposal be adapted during sampling} \item{proposalScale}{scaling factor for the proposal generation. Can be adapted if there is too much / too little rejection} + +\item{x}{Parameter to generate the exponential sequence for building intermediary distributions. Default value from Jeremiah et al. (2012)} + +\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"} + +\item{ess.limit}{Threshold value of effective sample size below which resampling is done. By default, the value is set to half the number of particles. To resample at each step, use a value >= the number of particles.} + +\item{ess.factor}{} + +\item{lastResample}{Iteration (starting from the end) at which particle resampling is forced. To deactivate this, set to a value < 0} + +\item{pars.lower}{Optional: vector contaiing the mimimum values of calibration parameters. Used to initialize the proposal function} + +\item{pars.upper}{Optional: vector contaiing the maximum values of calibration parameters. Used to initialize the proposal function} + +\item{mutate.method}{MCMC method for the resampling. Either "Metropolis" or "DE"} + +\item{b}{parameter for snooker update in DE} + +\item{diagnostics}{an optional function with diagnostics that are calculated on the particles during each SMC iteration} } \description{ Sequential Monte Carlo Sampler @@ -46,97 +86,17 @@ The SMC currently assumes that the initial particle is sampled from the prior. I # 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), +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) - -} - - - +SMCout <- smcSampler(bayesianSetup) +marginalPlot(SMCout, type = "d", singlePanel = T) } \author{ -Florian Hartig +Florian Hartig, Matthias Speich } From b87c0b6385ed85cd76e56d08a846bcaa222c9bb3 Mon Sep 17 00:00:00 2001 From: florianhartig Date: Tue, 20 Apr 2021 17:28:20 +0200 Subject: [PATCH 36/42] integration in runMCMC and help updates --- BayesianTools/R/SMC.R | 6 +- BayesianTools/R/mcmcRun.R | 67 +++++++------------ BayesianTools/inst/examples/SMCHelp.R | 17 +++-- .../man/getMetropolisDefaultSettings.Rd | 15 ----- BayesianTools/man/smcSampler.Rd | 17 +++-- .../tests/testthat/test-settingsDefault.R | 2 + 6 files changed, 52 insertions(+), 72 deletions(-) delete mode 100644 BayesianTools/man/getMetropolisDefaultSettings.Rd diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index 7d48edf..9f927d8 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -41,12 +41,12 @@ smcSampler <- function(bayesianSetup, m=7E-08, sampling="multinomial", ess.limit=NULL, - ess.factor = 0.95, # TODO document + ess.factor = 0.95, lastResample = 1, pars.lower=NULL, pars.upper=NULL, - mutate.method ="Metropolis", # TODO document - b=1e-04, # TODO document + mutate.method ="Metropolis", + b=1e-04, diagnostics = NULL, reference=NULL){ diff --git a/BayesianTools/R/mcmcRun.R b/BayesianTools/R/mcmcRun.R index 48d8e83..d37fec7 100644 --- a/BayesianTools/R/mcmcRun.R +++ b/BayesianTools/R/mcmcRun.R @@ -280,17 +280,17 @@ 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") if (settings$sampler == "AM") { - defaultSettings <- getMetropolisDefaultSettings() + defaultSettings <- applySettingsDefault(sampler = "Metropolis") defaultSettings$adapt <- TRUE } if (settings$sampler == "DR") { - defaultSettings <- getMetropolisDefaultSettings() + defaultSettings <- applySettingsDefault(sampler = "Metropolis") defaultSettings$DRlevels <- 2 } if (settings$sampler == "DRAM") { - defaultSettings <- getMetropolisDefaultSettings() + defaultSettings <- applySettingsDefault(sampler = "Metropolis") defaultSettings$adapt <- TRUE defaultSettings$DRlevels <- 2 } @@ -317,8 +317,6 @@ applySettingsDefault<-function(settings=NULL, sampler = "DEzs", check = FALSE){ message = TRUE) } - - if (settings$sampler == "DE"){ defaultSettings = list(startValue = NULL, iterations = 10000, @@ -408,13 +406,27 @@ applySettingsDefault<-function(settings=NULL, sampler = "DEzs", check = FALSE){ } 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 ) } @@ -556,33 +568,4 @@ getPossibleSamplerTypes <- function(){ ) return(out) -} - -#' Returns Metropolis default settings -#' @author Tankred Ott -#' @keywords internal -getMetropolisDefaultSettings <- function () { - 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 - ) - return(defaultSettings) -} - +} \ No newline at end of file diff --git a/BayesianTools/inst/examples/SMCHelp.R b/BayesianTools/inst/examples/SMCHelp.R index fcf4dc2..5704dc3 100644 --- a/BayesianTools/inst/examples/SMCHelp.R +++ b/BayesianTools/inst/examples/SMCHelp.R @@ -1,14 +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 - - +# creating a simple test BayesianSetup bayesianSetup <- createBayesianSetup(likelihood = generateTestDensityMultiNormal(), lower = rep(-10, 3), upper = rep(10, 3), parallel = "external") -SMCout <- smcSampler(bayesianSetup) +# 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/man/getMetropolisDefaultSettings.Rd b/BayesianTools/man/getMetropolisDefaultSettings.Rd deleted file mode 100644 index a2dbfe3..0000000 --- a/BayesianTools/man/getMetropolisDefaultSettings.Rd +++ /dev/null @@ -1,15 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/mcmcRun.R -\name{getMetropolisDefaultSettings} -\alias{getMetropolisDefaultSettings} -\title{Returns Metropolis default settings} -\usage{ -getMetropolisDefaultSettings() -} -\description{ -Returns Metropolis default settings -} -\author{ -Tankred Ott -} -\keyword{internal} diff --git a/BayesianTools/man/smcSampler.Rd b/BayesianTools/man/smcSampler.Rd index eb2b4f3..2328d47 100644 --- a/BayesianTools/man/smcSampler.Rd +++ b/BayesianTools/man/smcSampler.Rd @@ -82,19 +82,24 @@ The sampler can be used for rejection sampling as well as for sequential Monte C 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 - - +# creating a simple test BayesianSetup bayesianSetup <- createBayesianSetup(likelihood = generateTestDensityMultiNormal(), lower = rep(-10, 3), upper = rep(10, 3), parallel = "external") -SMCout <- smcSampler(bayesianSetup) +# 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{ diff --git a/BayesianTools/tests/testthat/test-settingsDefault.R b/BayesianTools/tests/testthat/test-settingsDefault.R index a98dd15..4aebab8 100644 --- a/BayesianTools/tests/testthat/test-settingsDefault.R +++ b/BayesianTools/tests/testthat/test-settingsDefault.R @@ -11,6 +11,8 @@ test_that("Default works in principle",{ settings <- list(iterations = 20000, adapt = T, DRlevels = 2, optimize = T, burnin=1000, adaptationInterval=10) applySettingsDefault(settings = settings, sampler = "Metropolis", check = FALSE) applySettingsDefault(settings = settings, sampler = "Metropolis", check = TRUE) + + applySettingsDefault(settings = settings, sampler = "DRAM", check = TRUE) } ) From 174a2026a6b7a4bb84096f8647efa58f9f34f5a9 Mon Sep 17 00:00:00 2001 From: mspeich Date: Wed, 21 Apr 2021 21:44:39 +0200 Subject: [PATCH 37/42] Expand parameter descriptions for SMC function --- BayesianTools/R/SMC.R | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index 9f927d8..1d3976c 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -5,22 +5,22 @@ #' @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 if new particles should be created at each iteration -#' @param resamplingSteps how many resampling (MCMC) steps between the SMC iterations -#' @param lastMutateSteps how many resampling (MCMC) steps after the SMC iterations +#' @param resamplingSteps how many resampling (MCMC) steps between the SMC iterations. A value of 20 or higher is recommended. +#' @param lastMutateSteps how many resampling (MCMC) steps after the SMC iterations to increase sample diversity. #' @param proposal optional proposal class #' @param exponents series of exponents to build the intermediate distributions #' @param adaptive 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 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 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 of effective sample size below which resampling is done. By default, the value is set to half the number of particles. To resample at each step, use a value >= the number of particles. +#' @param ess.limit Threshold value of effective sample size below which resampling is done (fraction of effective sample size). By default, the value is set to half the number of particles. To resample at each step, use a value >= the number of particles. #' @param ess.factor #' @param lastResample Iteration (starting from the end) at which particle resampling is forced. To deactivate this, set to a value < 0 -#' @param pars.lower Optional: vector contaiing the mimimum values of calibration parameters. Used to initialize the proposal function -#' @param pars.upper Optional: vector contaiing the maximum values of calibration parameters. Used to initialize the proposal function -#' @param mutate.method MCMC method for the resampling. Either "Metropolis" or "DE" -#' @param b parameter for snooker update in DE +#' @param pars.lower Optional: vector containing the mimimum values of calibration parameters. Used to initialize the proposal function +#' @param pars.upper Optional: vector containing the maximum values of calibration parameters. Used to initialize the proposal function +#' @param mutate.method MCMC method for the resampling. Either "Metropolis" (default) or "DE" +#' @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. From 84e1ead6076a5ecbf9895ca9834532d3e322ad3b Mon Sep 17 00:00:00 2001 From: mspeich Date: Wed, 21 Apr 2021 21:51:26 +0200 Subject: [PATCH 38/42] Add script to plot SMC output Add the script used to generate Fig. 3 in Speich et al. --- .../7_plot_results/process_smc_output.r | 317 ++++++++++++++++++ 1 file changed, 317 insertions(+) create mode 100644 Publications/SpeichEtAl-SMCForModelCalibration/7_plot_results/process_smc_output.r 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 From 4bbfb4eed52f77ca39bb9d3e6baa5cb6d10900b5 Mon Sep 17 00:00:00 2001 From: florianhartig Date: Thu, 22 Apr 2021 22:30:06 +0200 Subject: [PATCH 39/42] Update SMC.R --- BayesianTools/R/SMC.R | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index 1d3976c..c73e48a 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -4,7 +4,7 @@ #' @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 if new particles should be created at each iteration +#' @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 lastMutateSteps how many resampling (MCMC) steps after the SMC iterations to increase sample diversity. #' @param proposal optional proposal class @@ -30,7 +30,7 @@ smcSampler <- function(bayesianSetup, initialParticles = 1000, iterations = 10, - resampling = T, + resampling = NULL, resamplingSteps = 2, lastMutateSteps = 5, proposal = NULL, @@ -52,6 +52,8 @@ smcSampler <- function(bayesianSetup, # 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() @@ -239,7 +241,6 @@ smcSampler <- function(bayesianSetup, ess <- inter.out$ess doResample <- inter.out$doResample - info$ess.vec[icount] <- ess info$exponents[icount] <- curExp @@ -247,7 +248,7 @@ smcSampler <- function(bayesianSetup, curExp <- 1 } - # Re?sampling step + # 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 @@ -273,9 +274,6 @@ smcSampler <- function(bayesianSetup, # 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)))) @@ -294,7 +292,7 @@ smcSampler <- function(bayesianSetup, 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(curExp, nrow(unique(particles)))) @@ -327,8 +325,6 @@ smcSampler <- function(bayesianSetup, # info$resamplingAcceptance[(icount),] <- mutate.out$acceptance # } - - oldweights <- weights info$ess.vec[icount] <- ess From 63f42ffa3035e77fcb4f1a9a7c2cfea0b59210d4 Mon Sep 17 00:00:00 2001 From: florianhartig Date: Thu, 22 Apr 2021 22:46:01 +0200 Subject: [PATCH 40/42] updated parameter sequence in help --- BayesianTools/R/SMC.R | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index c73e48a..168b4a5 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -4,25 +4,27 @@ #' @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 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 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 exponents series of exponents to build the intermediate distributions -#' @param adaptive 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. 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 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 of effective sample size below which resampling is done (fraction of effective sample size). By default, the value is set to half the number of particles. To resample at each step, use a value >= the number of particles. #' @param ess.factor #' @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 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 mimimum values of calibration parameters. Used to initialize the proposal function #' @param pars.upper Optional: vector containing the maximum values of calibration parameters. Used to initialize the proposal function #' @param mutate.method MCMC method for the resampling. Either "Metropolis" (default) or "DE" #' @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. #' @export From ccca96f1793f47a140e38ae57c618d5b7bf873fc Mon Sep 17 00:00:00 2001 From: florianhartig Date: Sat, 10 Jul 2021 09:51:52 +0200 Subject: [PATCH 41/42] Update SMC.R --- BayesianTools/R/SMC.R | 374 +++++++++++++++++++++--------------------- 1 file changed, 187 insertions(+), 187 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index 168b4a5..0c16891 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -1,10 +1,10 @@ #' SMC sampler #' @author Florian Hartig, Matthias Speich #' @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. +#' @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 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 @@ -12,56 +12,56 @@ #' @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 of effective sample size below which resampling is done (fraction of effective sample size). By default, the value is set to half the number of particles. To resample at each step, use a value >= the number of particles. -#' @param ess.factor +#' @param ess.factor #' @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 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 mimimum values of calibration parameters. Used to initialize the proposal function -#' @param pars.upper Optional: vector containing the maximum values of calibration parameters. Used to initialize the proposal function -#' @param mutate.method MCMC method for the resampling. Either "Metropolis" (default) or "DE" +#' @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. +#' @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, +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, + 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=NULL, - ess.factor = 0.95, + ess.factor = 0.95, lastResample = 1, - pars.lower=NULL, - pars.upper=NULL, - mutate.method ="Metropolis", - b=1e-04, + 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") - + + # 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() # 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). @@ -73,10 +73,10 @@ smcSampler <- function(bayesianSetup, setup <- checkBayesianSetup(bayesianSetup) - + ### InitialParticles #TODO documentation - + if(class(initialParticles)[1] == "numeric"){ initialParticles = bayesianSetup$prior$sampler(initialParticles) importanceDensity = bayesianSetup$prior$density @@ -94,32 +94,32 @@ smcSampler <- function(bayesianSetup, initialParticles <- initialParticles$particles } } - + if (any(is.infinite(setup$prior$density(initialParticles)))) stop("initialParticles outside prior range") - + particles <- initialParticles rejectionRate = 0 particleSize = nrow(initialParticles) if(is.null(ess.limit)){ess.limit <- round(particleSize) * 0.5} - + weights <- oldweights <- oldInter <- rep(0, particleSize) - + posterior = matrix(nrow = particleSize, ncol = 3) numPar <- ncol(initialParticles) - + if(mutate.method %in% c("Metropolis", "adaptive")){ if (is.null(proposal)) proposalGenerator = createProposalGenerator(rep(40,numPar)) } - + usedUp = 0 - + #### TEST #proposalScale <- runif(nrow(particles), min = 1e-06, max = 1) Z <- particles - - # WEIGHT SETUP - + + # 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)){ @@ -128,7 +128,7 @@ smcSampler <- function(bayesianSetup, # 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 @@ -136,41 +136,41 @@ smcSampler <- function(bayesianSetup, 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) + 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") @@ -179,20 +179,20 @@ smcSampler <- function(bayesianSetup, 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) @@ -200,41 +200,41 @@ smcSampler <- function(bayesianSetup, max.pars <- pmax(max.pars, cur.max) if (numPar == 1) particles = matrix(particles, ncol = 1) - + # 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 - + # 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 - + + # 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 - + #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-1), posteriorValues = posteriorValues, importanceValues = importanceValues, oldInter = oldInter, curWeights = weights, curExp = curExp) curExp <- inter.out$newExp @@ -242,42 +242,42 @@ smcSampler <- function(bayesianSetup, 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 + + # 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 | 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 + # Scale particles #print(c("nrow unique particles", nrow(unique(particles)))) #print(c("curExp", curExp)) #if(nrow(unique(particles)) < 1000){ @@ -285,7 +285,7 @@ smcSampler <- function(bayesianSetup, #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 @@ -295,14 +295,14 @@ smcSampler <- function(bayesianSetup, if(length(proposalScale) > 1){ proposalScale <- mutate.out$proposalScale } - + #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 @@ -310,14 +310,14 @@ smcSampler <- function(bayesianSetup, # 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) @@ -326,26 +326,26 @@ smcSampler <- function(bayesianSetup, # importanceValues <- mutate.out$importanceValues # info$resamplingAcceptance[(icount),] <- mutate.out$acceptance # } - + oldweights <- weights info$ess.vec[icount] <- ess - - icount <- icount + 1 + + 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 @@ -355,19 +355,19 @@ smcSampler <- function(bayesianSetup, 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. @@ -379,9 +379,9 @@ smcSampler <- function(bayesianSetup, 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, settings = settings, @@ -389,15 +389,15 @@ smcSampler <- function(bayesianSetup, posterior = posterior, info = info ) - + class(out) <- c("smcSampler", "bayesianOutput") return(out) - + } #' Residual Resampling -#' +#' #' @keywords internal residualResampling <- function(weights){ weights <- exp(weights) @@ -405,12 +405,12 @@ residualResampling <- function(weights){ 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) @@ -424,7 +424,7 @@ systematicResampling <- function(weights){ 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) @@ -432,24 +432,24 @@ systematicResampling <- function(weights){ 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 -#' -#' +#' 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"){ @@ -463,14 +463,14 @@ resample <- function(weights, method = "multinomial"){ } -#' Auxiliary function for adaptive SMC algorithm -#' -#' +#' 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 @@ -478,16 +478,16 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI 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") @@ -510,7 +510,7 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI 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 @@ -522,7 +522,7 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI } - + if(any(is.infinite(tryWeights)) | is.infinite(try.ess)){ doResample <- TRUE tryExp <- curExp @@ -532,51 +532,51 @@ beta.search <- function(ess, target.ess, posteriorValues, importanceValues, oldI } 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 @@ -584,7 +584,7 @@ mutate <- function(setup, particles, proposalGenerator, posteriorValues, importa 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)){ @@ -597,26 +597,26 @@ mutate <- function(setup, particles, proposalGenerator, posteriorValues, importa 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{ + + } 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)){ @@ -629,9 +629,9 @@ mutate <- function(setup, particles, proposalGenerator, posteriorValues, importa } numPar <- ncol(particles) randVector <- runif(numPar,-b,b) - particlesProposals[part,] <- particles[part,] + (particleDiff * scale.factors[part]) + randVector + particlesProposals[part,] <- particles[part,] + (particleDiff * scale.factors[part]) + randVector r_extra[part] <- 0 - + # if(part %in% 1:6){ # print(c("Part ", part)) # print("Particle") @@ -643,23 +643,23 @@ mutate <- function(setup, particles, proposalGenerator, posteriorValues, importa # print(c("Part NO snooker", part, r_extra[part])) #} } - } + } } - - proposalPosteriors <- setup$posterior$density(particlesProposals) + + 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. @@ -667,7 +667,7 @@ mutate <- function(setup, particles, proposalGenerator, posteriorValues, importa 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") @@ -676,29 +676,29 @@ mutate <- function(setup, particles, proposalGenerator, posteriorValues, importa print(posteriorValues) print("r_extra") print(r_extra) - + } - - particles[accepted, ] = particlesProposals[accepted, ] + + 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 @@ -712,18 +712,18 @@ mutate <- function(setup, particles, proposalGenerator, posteriorValues, importa #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))) @@ -733,20 +733,20 @@ mutate <- function(setup, particles, proposalGenerator, posteriorValues, importa } 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){ @@ -754,13 +754,13 @@ mutate <- function(setup, particles, proposalGenerator, posteriorValues, importa } 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) @@ -769,7 +769,7 @@ scale.particles <- function(particle.row, min.pars, max.pars){ #' Calculate rao diversity -#' +#' #' @keywords internal rao.div <- function(scaled.particles){ # Determine unique particles, their sum and pairwise difference @@ -784,18 +784,18 @@ rao.div <- function(scaled.particles){ } }, 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) } @@ -805,5 +805,5 @@ rao.div <- function(scaled.particles){ # largeRatio <- function(a,b){ # # Determine the larger of the two numbers -# -# } \ No newline at end of file +# +# } From 65d712de657e77c40564db292dc940583173f20b Mon Sep 17 00:00:00 2001 From: Florian Hartig Date: Sun, 8 Mar 2026 15:41:05 +0100 Subject: [PATCH 42/42] updates --- BayesianTools/R/SMC.R | 30 +++-- BayesianTools/man/smcSampler.Rd | 149 +++++++++--------------- BayesianTools/tests/testthat/test-SMC.R | 2 +- 3 files changed, 77 insertions(+), 104 deletions(-) diff --git a/BayesianTools/R/SMC.R b/BayesianTools/R/SMC.R index cb97027..4947125 100644 --- a/BayesianTools/R/SMC.R +++ b/BayesianTools/R/SMC.R @@ -1,19 +1,18 @@ #' SMC sampler -#' @author Florian Hartig, Matthias Speich -#' @description Sequential Monte Carlo Sampler +#' +#' @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 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 of effective sample size below which resampling is done (fraction of effective sample size). By default, the value is set to half the number of particles. To resample at each step, use a value >= the number of particles. -#' @param ess.factor +#' @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 @@ -28,6 +27,9 @@ #' #' @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 +#' #' @export #' @example /inst/examples/SMCHelp.R smcSampler <- function(bayesianSetup, @@ -43,7 +45,7 @@ smcSampler <- function(bayesianSetup, x=3.11, m=7E-08, sampling="multinomial", - ess.limit=NULL, + ess.limit=0.5, ess.factor = 0.95, lastResample = 1, pars.lower=NULL, @@ -72,7 +74,7 @@ smcSampler <- function(bayesianSetup, 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) @@ -102,7 +104,7 @@ smcSampler <- function(bayesianSetup, particles <- initialParticles rejectionRate = 0 particleSize = nrow(initialParticles) - if(is.null(ess.limit)){ess.limit <- round(particleSize) * 0.5} + ess.limit.abs = round(particleSize) * ess.limit weights <- oldweights <- oldInter <- rep(0, particleSize) @@ -237,8 +239,16 @@ smcSampler <- function(bayesianSetup, 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-1), 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 @@ -258,7 +268,7 @@ smcSampler <- function(bayesianSetup, # 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 | icount == (length(exponents) - lastResample) | doResample){ + if(ess < ess.limit.abs | icount == (length(exponents) - lastResample) | doResample){ oldExp <- curExp 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 index 77ff58c..369cc14 100644 --- a/BayesianTools/tests/testthat/test-SMC.R +++ b/BayesianTools/tests/testthat/test-SMC.R @@ -19,7 +19,7 @@ test_that("resampling works", { weights = log(c(1,2,3)/6) size = 10000 - x1 = replicate(size, resample(weights)) + x1 = replicate(size, BayesianTools:::resample(weights)) dif = table(x1) - exp(weights) * size * 3 expect_lt(sum(abs(dif)), 300)