R/fpdust.R

In rhobis/FPdust: An Implementation of the FPDUST and PPS FPDUST Spatial Sampling Designs for Finite Populations

#' Select a sample from a finite population using FPDUST procedure.
#'
#' Select a sample of size n from a finite population given a spatial
#' contingency matrix, the sample size \code{n}, and a parameter \code{beta}.
#'
#' @param spMat a matrix representing spatial contingency among population units,
#'         with dimension \code{NxN} and whose values are 1's (adjacent units) or 0's (otherwisex)
#' @param n integer scalar representing the sample size
#' @param beta a real number <1 representing a modifier of the probability of
#'             selection for units that are adjacent to already selected units
#'
#' @return a logical vector indicating which units have been selected
#'         (TRUE if selected, FALSE otherwise)
#'
#'
#' @examples
#' # List of contiguous units ---
#' cl <- list( c(2, 3), c(1,4), c(1,6), c(2,6), NULL, c(3,4))
#' # Generate spatial contingency matrix ---
#' cm <- spatialContingencyMatrix(cl)
#'
#' ### Draw a FPDUST sample ---
#' fpdust(cm, n= 2, beta=0.5)
#'
#'
#'
#' @export


fpdust <- function(spMat,n,beta){

    ## Check input
    if(!is.matrix(spMat)) stop("Argument spMat should be anobject of class matrix!")
    if(ncol(spMat) != nrow(spMat)) stop("spMat should be a square matrix!")
    if( !is.atomic(n)) stop("Argument n must be a numeric scalar")
    if( !is.numeric(n) ) stop("Argument n must be numeric!")
    if( n<2) stop("The sample size n should be >=2")
    if( beta>1) stop("The parameter beta should be <1")


    if(beta >= 1) stop('The argument beta should be strictly less than 1')
    n <- floor(n)
    N <- nrow(spMat)

    if(N < n) stop('n is greater than population size')
    diag(spMat) <- 0 #to be sure units are not considered adjacent to themselves
    sv <- rep(FALSE,N) #vector of sampled units
    ct <- rep(FALSE,N) #vector of adjacent units

    sv[sample(1:N, 1)] <- TRUE
    ct <- as.logical(spMat[sv,])

    for(k in 2:n){
        p <- (1-sv) * (1 - beta*ct) / (N - k + 1 - sum(ct)*beta)
        # print(sum(p))
        i <- sample(1:N,1,prob=p)
        sv[i] <- TRUE
        ct <- ct | spMat[i,]
        ct <- ct * xor(ct,sv)
    }

    return(sv)

}