weierstrass: Combine independent posterior samples from multiple subsets
In wwrechard/weierstrass: Combining posterior samples from multiple subsets.

Description Usage Arguments Details Value See Also Examples

View source: R/weierstrass.r

Combine independent posterior samples from multiple subsets

weierstrass(Samples, num.sets = NULL, para.dim = NULL, method = "reject",
  kernel = "uniform", kernel.corr = TRUE, accept = 0.1, weight = TRUE,
  average = TRUE, matching = TRUE, resampling = TRUE,
  resampling.corr = FALSE, display = TRUE, level = 1, sets = 1)

`Samples`	A list containing posterior samples from each subset. Samples[[i]] contains a matrix of posterior samples for i-th set. Each line of the matrix is a sample and each column corresponds to a parameter.
`num.sets`	The number of subsets. If not specified, the length of Samples will be used.
`para.dim`	The number of parameters. If not specified, The number of columns of Samples[[1]] will be used.
`method`	Methods for combining the subset samples. Default is "reject" for rejection sampling. Alternative could be "importance" for importance sampling.
`kernel`	The kernel function used for rejection sampling and importance sampling. The default is "uniform". Alternative could be "gaussian". Importance sampling always uses "gaussian" no matter what value is specified for kernel.
`kernel.corr`	Indicator of whether a correlated kernel should be used for each subset. If TRUE, the algorithm will use a correlated kernels where the correlation equal to the sample correlation of posteriors on each subset. If FALSE, independent kernels are used. Default is TRUE.
`accept`	The acceptance rate for rejection sampling and importance resampling. For more information please refer to the details. Default value is 0.1.
`weight`	Indicator of whether subsets should be weighted. If TRUE, inverse variance weighting will be used. Depending on kernel.corr, the weights could be either matrices or scalars. Default is TRUE.
`average`	Indicator of how samples are combined for initial processing. If TRUE, the subset samples will be (weighted) averaged before entering the rejection stage. If FALSE, mixture will be used instead, i.e., randomly (or weighted) select one sample among all subsets as the candidate for rejection sampling. The parameter is always set to TRUE for importance sampling. For more information please refer to the details. Default is TRUE.
`matching`	Number of samples regenerated at each iteration (and final output). If TRUE, the algorithm will match the number of samples generated at each iteration to be the same as the number of input posterior samples (on one subset). If a number is given, the algorithm will use that number for each iteration. If FALSE, the algorithm will adjust its acceptance rate at each iteration to ensure the final number of combined samples equal to the input size times acceptance rate. For more information please refer to the details. Default is TRUE.
`resampling`	Indicate whether resampling should be used for importance sampling. Default is TRUE.
`resampling.corr`	If resampling is set to be TRUE, resampling.corr determines whether a joint kernel should be used for resampling, otherwise independent kernels are used. Default is FALSE.
`display`	Indicate whether processing information should be printed. Default is TRUE.

Weierstrass is the implementation of the Weierstrass rejection sampler proposed in the paper. Weierstrass sampler is a 'divide-conquer-combine' type parallel sampler, including three distinct samplers. Different from the refining sampler and the sequential rejection sampler, Weierstrass rejection sampler is a pure post-processing algorithm in the sense that it directly works on posterior samples obtained from multiple subsets, while the other two also relies on the subset sampling procedure and the likelihood. The algorithm applies rejection sampling and importance sampling on the formula described in Theorem 3 in the paper.

To combine the subset posterior samples, the algorithm adopts the 'pairwise-combining' strategy, i.e, it will first combine the subset pairwisely to obtain half numbers of new subsets and then repeat the procedure until obtaining the final one. Such procedure, though, adding a logarithm factor to the original run time, achieves substantially better result. Based on this process, the parameter "matching" is introduced as for determining how many samples should be preserved (or regenerated) at each iteration. If "matching" is set to be FALSE, the algorithm will increase the acceptance rate at each iteration in order to obtain enough numbers of combined samples at the final output. If "matching" is set to be TRUE or any number, the algorithm will use the same acceptance rate at each iteration but regenerate samples to match the number required. For the final output, the number of combined samples will be equal to "matching" (if a number is assigned to "matching") or the number of input samples (If "matching" is equal to TRUE).

Subset samples are required to be initially combined before entering the rejection sampling phase. Two different ways for initial combining are introduced. The first is (weighted) average(if "average" is TRUE, default) and the second is mixture (otherwise). Typically, weighted average achieves better accuracy, which is thereby recommended for most cases. However, there are exceptions. First, when the support of posterior samples are discontinuous. For example, integers or confined regions, in which case, any averaging procedure will fail. Mixture is a good alternative for these examples as it will always be in the correct region. Second, when the posterior is multimodal and each parameter is likely to stuck at one mode. For example, Gaussian mixture model without label switching. In this case, using mixture is likely to recover the multimodel shape of posteriors at the cost of losing (a certain level) accuracy. It is worth noting that when the parameter "accept" is set to be 1, the algorithm will combine the subset posteriors just by (weighted) averaging or mixture, which might serve for other interest.

"accept" is the only parameter that needs to be manually determined by the user. It determines the acceptance rate for rejection sampling (and importance sampling). The bandwidth will also be tuned automatically by the algorithm based on the value of acceptance rate. It is crucial taht "accept" have impacts on both the accuracy and the complexity of the algorithm. When "matching" is set to be TRUE, the run time of the algorithm will be proportional to 1/accept. When the "matching" is set to be FALSE, the number of actual obtained combined posterior samples is equal to the number of input posterior samples times "accept". In both cases, the value of "accept" cannot be set too small. In addition, even in the scenario when "matching" is TRUE, "accept" will also subtly affect the effective sample size. Therefoe, we recommend the value of "accept" to be chosen within (0.01, 0.1) or just set equal to 2/(number of subsets). The default is 0.1.

A matrix contain combined posterior samples. Each row is a sample and each column corresponds to a parameter.

logitTest, BinTest

## Not run: Two examples are provided as separate functions. \code{\link{logitTest}} is an example for logistic regression which depends on the package "BayesLogit" and \code{\link{BinTest}} is an example for binomial distributions.
## Not run: A simple example of linear regression is given as follows.
library('MASS')
p = 5
n = 10000
m = 20
V = n/m
beta = matrix(c(1,-1,2,4,0),p,1)
X = matrix(rnorm(n*p), n, p)
err = matrix(rnorm(n),n,1)
Y = X%*%beta + err
piror.mu = matrix(rep(0,p),p,1)
prior.var = diag(rep(10,p),p,p)

Samples = list()
for(i in 1:m){
Samples[[i]] = mvrnorm(n = 10000, mu = solve(t(X[(V*(i-1)+1):(V*i),])%*%X[(V*(i-1)+1):(V*i),] + solve(prior.var)/m)%*%t(X[(V*(i-1)+1):(V*i),])%*%Y[(V*(i-1)+1):(V*i),], Sigma = solve(t(X[(V*(i-1)+1):(V*i),])%*%X[(V*(i-1)+1):(V*i),] + solve(prior.var)/m))
}
true.posterior = mvrnorm(n = 10000, mu = solve(t(X)%*%X + solve(prior.var)/m)%*%t(X)%*%Y, Sigma = solve(t(X)%*%X + solve(prior.var)/m))
CombSample = weierstrass(Samples)

par(mfrow = c(2,3))
for(i in 1:5){
   plot(density(true.posterior[,i]),col = 'red',main = paste('posterior for beta',toString(i)))
   lines(density(CombSample[,i]),col = 'blue')
   legend('topleft', c('True posterior', 'Weierstrass'), col = c('red','blue'), lty = c(1,1))
}