R/BinTest.r

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

#' A test function provide binomial data for testing Weierstrass rejection sampling (and comparing to the fullset posterior/averaging/weighted averaging combiners).
#'
#' @param prob The parameter for binomial distribution. Default is 0.001.
#' @param samp The posterior sample size on each subset. Default is 20000.
#' @param n The total sample size. Default is 10000.
#' @param m Number of subsets. Default is 20.
#' @param accept The acceptance rate for the Weierstrass rejection sampling. Default is 0.2.
#' @param draw Indicate whether the result should be plotted.
#'
#' @return A list containing several components.
#' \enumerate{
#' \item posteriors: a list containing all subset posterior samples. Input for weierstrass rejection samping.
#' \item true.posterior: a matrix containing posterior samples drawn with full data set.
#' \item CombSample.weight: a matrix containing combined samples generated by the Weierstrass rejection sampler
#' \item CombSample.unweight: a matrix containing combined samples generated by the unweighted Weierstrass rejection sampler
#' \item weight.ave: a matrix containing combined samples via inverse-variance weighted averaging.
#' \item ave: a matrix containing combined samples via simple averaging.
#' }
#'
#' @details The function generates data from a binomial distribution and make use of Weierstrass rejection sampler, weighted average and averge to combine the subset posterior samples.Both the weighted and unweighted weierstrass sampler will be used, and the results are saved for further processing.
#'
#' @seealso \code{\link{weierstrass}} for the details of weierstrass rejection sampling. \code{\link{logitTest}} for another test on the logistic regression.
#' 
#' @examples
#' \dontrun{BinTest()}
#' \dontrun{BinTest(prob = 0.49)}
#' \dontrun{BinTest(prob = 0.99)}
#' 
BinTest <- function(prob = 0.001, samp = 20000, n = 20000, m = 20, accept = 0.1, draw = TRUE){
    p = 2
    V = floor(n/m)
    N = rep(V,m)
    N = c(0,N)
    N = cumsum(N)
    
    prior = rep(0.01,p)
    sub.prior = prior/m

    pt = c(prob, 1-prob)
  
    count = rep(0,p)
    pt.cum = cumsum(pt)
    experiment = runif(n)

    count = matrix(0,m,p)
    for(j in 2:(m+1)){
        expe = experiment[(N[j-1]+1):N[j]]
        for(i in 1:p){
            count.temp = expe<=pt.cum[i]
            expe = expe[!count.temp]
            count[j-1,i] = sum(count.temp)
        }
        count[j-1,] = count[j-1,]+sub.prior
    }

    posterior<-list()
    for(j in 1:m){
        posterior[[j]] = matrix(rbeta(samp,sub.prior[1]+count[j,1]+1, sub.prior[2]+sum(count[j,-1])+1),samp,p-1)
    }
    true.posterior = matrix(rbeta(samp,prior[1]+apply(count,2,sum)[1]+1, prior[2]+sum(apply(count,2,sum)[-1])+1),samp,p-1)

    #combine via (weighted) averaging
    ave1 = ave2 = matrix(0, samp, 1)
    icov = list()
    bigc = matrix(0,1,1)
    for(i in 1:m){
        icov[[i]] = solve(cov(posterior[[i]]))
        bigc = bigc + icov[[i]]
    }
    for(i in 1:m){
        ave1 = ave1 + posterior[[i]]%*%icov[[i]]%*%solve(bigc)
        ave2 = ave2 + posterior[[i]]
    }
    ave2 = ave2/m

    #combine via Weierstrass rejection sampler
    cat(paste('starting the weierstrass rejection sampling, acceptance rate is',toString(accept),'\n\n'))
    CombSample1 = weierstrass(posterior, accept = accept)
    
    cat('completed.\n\n')
    cat(paste('starting the unweighted weierstrass rejection sampling, acceptance rate is',toString(accept),'\n\n'))
    CombSample2 = weierstrass(posterior, weight = FALSE, accept = accept)
    cat('completed.\n\n')
    
    #plot the result
    if(draw == TRUE){
        rows = floor(sqrt(p-1))
        cols = ceiling((p-1)/rows)
        par(mfrow = c(rows, cols))
        for(i in 1:(p-1)){
            plot(density(true.posterior[,i]),col='red',main = paste('Posterior for the parameter'))
            lines(density(CombSample1[,i]),col = 'blue')
            lines(density(CombSample2[,i]),col = 'yellow',lty = 2)
            lines(density(ave1[,i]),col = 'green')
            lines(density(ave2[,i]),col = 'pink')
            abline(v = prob, col = 'black', lty = 2)
            if(prob>0.5) posi = 'topleft'
            else posi = 'topright'
            legend(posi,c('Fullset','Weierstrass', 'Unweight weiers' ,'Weighted ave','average','true value'),col = c('red','blue','yellow','green','pink','black'),lty = c(1,1,2,1,1,1))
        }
    }
    return(list(Samples = posterior, true.posterior = true.posterior, CombSample.weight = CombSample1, CombSample.unweight = CombSample2, Weight.ave = ave1, ave = ave2))
}