R/simulations.R

In DistributionalDDA: Distributional Data Analysis in Bayes Spaces

# This files contains simulation functions generating densities

#simulate = function(object, nsim = 1, seed = NULL, ...) UseMethod("simulate",x)


# ## uniform reference measure
# ReferenceMeasureUniform = function(a,b,nsim=10000,sample=runif(nsim,a,b)) {
#   structure(list(a=a,b=b,sample=sample),
#             class=c("ReferenceMeasureUniform","ReferenceMeasure"))
# } ### IN BAYES_SPACES.R
#' simulation method for the uniform referencemeasure
#'
#' @param object the ReferenceMeasureUniform object
#' @param nsim the number of simulations
#' @param seed optional seed, if set it is used and reported as
#'     attribute in the simulation (use with care).
#' @param ... further unused generic arguments
#'
#' @returns a vector of length nsim with of uniformly distributed
#'     random variables.
#' @export
#' @importFrom stats simulate
#' @importFrom stats runif
#' @seealso [stats::simulate]
#' @examples
#' set.seed(4)
#' simulate(StdReferenceMeasureUnif01,10,seed=2)
#' simulate(StdReferenceMeasureUnif01,10)
#' simulate(StdReferenceMeasureUnif01,10,seed=2)
#' simulate(StdReferenceMeasureUnif01,10)
simulate.ReferenceMeasureUniform = function(object,nsim=1L,seed=NULL,...) {
  referenceMeasure = getReferenceMeasure(object)
  a <- referenceMeasure$a
  b <- referenceMeasure$b
  seededSimulation(runif(nsim,a,b),seed=seed)
}

# Helper function to assign a seed
seededSimulation = function(result,seed=NULL) {
    if( is.null(seed) ) {
        result
    } else {
    #    oState = .Random.seed
        set.seed(seed)
        tryCatch(structure(result,seed=seed) #,
                 # finally=.Random.seed <<- oState
                 # .Random.seed <<- oState # it is not allowed by CRAN policies to change .Random.seed
                 )
    }
}

#' Generates simulations from the Reference Measure seen as distribution.
#'
#' The simulations are done by applying the transformation to simulations
#' on the underlying reference measure
#'
#' @param object the ReferenceMeasureUniform object
#' @param nsim the number of simulations
#' @param seed optional seed, if set it is used and reported as
#'     attribute in the simulation (use with care).
#' @param ... further unused generic arguments
#'
#' @returns a vector of length nsim with of random variables with the
#'     given distribution
#'
#'
#' @export
#' @examples
#' set.seed(4)
#' simulate(StdReferenceMeasureN01,10,seed=2)
#' simulate(StdReferenceMeasureN01,10)
#' simulate(StdReferenceMeasureN01,10,seed=2)
#' simulate(StdReferenceMeasureN01,10)
#'
simulate.ReferenceMeasureTransformed=function(object, nsim = 1, seed = NULL, ...) {
    structure(
        object$qf(simulate(object$reference,nsim=nsim,seed=seed,...)),
        seed=seed)
}

#' Generates simulations from the Discrete Reference Measure seen as distribution.
#'
#' The simulation is done by the r function sample
#'
#' @param object the ReferenceMeasureDiscrete object
#' @param nsim the number of simulations
#' @param seed optional seed, if set it is used and reported as
#'     attribute in the simulation (use with care).
#' @param ... further unused generic arguments
#'
#' @returns a vector (or list) of length nsim with of random variables with the
#'     given distribution
#'
#'
#' @export
#' @examples
#' set.seed(4)
#' aDiscreteMeasure = ReferenceMeasureDiscrete(c(0.2,0.3,0.4),X=c("A","B","C"))
#' simulate(aDiscreteMeasure,10,seed=2)
#' simulate(aDiscreteMeasure,10)
#' simulate(aDiscreteMeasure,10,seed=2)
#' simulate(aDiscreteMeasure,10)
#'
simulate.ReferenceMeasureDiscrete=function(object, nsim = 1, seed = NULL, ...) {
    if( !is.null(seed) ) {
    #     oState = .Random.seed
        set.seed(seed)
    }
    if(is.function(object$sampler)) {
        x = object$sampler(nsim)
    } else {
        x = object$X[sample(1:length(object$X),nsim,replace=TRUE,prob=object$p)]
    }
    # if( !is.null(seed) ) {
      # .Random.seed <<- oState # it is not allowed by CRAN policies to change .Random.seed
    # }
    structure(x,seed=seed)
}

#' Simulates from a ReferenceMeasureCartesianProduct
#'
#' @param object the ReferenceMeasureCartesianProduct object
#' @param nsim the number of simulations
#' @param seed optional seed, if set it is used and reported as
#'     attribute in the simulation (use with care).
#' @param ... further unused generic arguments
#'
#' @returns a dataframe with nsim rows with simulations of the product measure.
#' @export
#' @seealso [stats::simulate]
#' @examples
#' set.seed(4)
#' CP = cartesianProduct(StdReferenceMeasureUnif01,StdReferenceMeasureN01)
#' simulate(CP,2,seed=17)
#' simulate(CP,2)
#' simulate(CP,2,seed=17)
#' simulate(CP,2)
simulate.ReferenceMeasureCartesianProduct <- function(object,nsim=1L,seed=NULL,...) {
    seededSimulation({
        as.data.frame(lapply(object,simulate,nsim=nsim,...))
    },seed=seed)
}



#' Simulate from the distribution described by the BData objects
#'
#' @note **Not yet implemented.** Calling this method currently stops with an
#'     error. It is registered so that `simulate()` on a `BData` object gives
#'     an informative message rather than a method-dispatch failure. The
#'     documented parameters and return value describe the planned interface.
#'
#' @param object a dataset of distributions from a BayesSpace
#' @param nsim number of simulations
#' @param seed the seed to be used. If seed is given the resulting
#'     object has a seed attribute, with the seed given (use with care).
#' @param ... futher generic parameters e.g. for the method of
#'     simulation used.
#' @param returnBDataSamples a boolean, if returnBDataSamples==FALSE
#'     'object' may only contain 1 distribution and the simulation of
#'     that distribution is returned directly. If TRUE the simulation
#'     is done for each of them and the results is returned as
#'     BDataSamples.
#'
#' @returns (Planned, once implemented:) A BDataSample object with each entry
#'     nsim simulations and one simulation for each distribution in the BData
#'     object. Currently the function always throws an error.
#' @export
#' @method simulate BData
#'
simulate.BData = function(object, nsim = 1L, seed = NULL, ...,returnBDataSamples=numberOfCases(object)>1) {
    stop("Simulation of BData not yet implemented")
}

# TODO
# #' @describeIn simulate.BData specific method for "BDataFrequencyHistograms"
# #' objects, as these contain information about sample size
# #' @export
# #' @method simulate BDataFrequencyHistograms
# simulate.BDataFrequencyHistograms = function(object, nsim=1L, seed = NULL, ..., returnBDataSamples=numberOfCases(object)>1) {
#   stop("Simulation of BDataFrequencyHistograms not yet implemented")
# }

#' A distribution of distributions based a distribution family with
#' random hyperparameters.
#'
#'
#' @param density A function that can be called with a vectorized x of
#'     data values of the BayesSpace and parameters, to calculate a
#'     density.
#' @param basis A BayesSpaceBasis (or an object that can be converted
#'     to one) of BayesSpace containing the values of the distribution
#'     (of distributions).
#' @param paramSimulator A function to be called without arguments to
#'     simulate one set of parameters.
#' @param ... parameters to be given to parametersimulator
#' @param .hyperparameters the list of parametrs to be given to the
#'     the parameterSimulator
#' @param simulator A function that can be called with an integer n
#'     and the parameters to create n random simulations.
#'
#' This function creates a distribution, which can be used to simulate
#' distributional objects according to the random hyperparameters.
#'
#' @returns A B2DistHyper object encoding the distribution.
#'
#' @export
#  @family B2DistHyper
#' @examples
#' D = B2DistHyper(dbeta,B2Unif01Basis5,function() rexp(2),simulator=rbeta)
#' simulate(D,4,seed=47)
#'
B2DistHyper = function(
                       density,
                       basis,
                       paramSimulator,
                       ...,
                       simulator=NULL,
                       .hyperparameters=list(...)) {
    basis = as.BayesSpaceBasis(basis)
    stopifnot(is.function(density))
    stopifnot(is.function(paramSimulator))
    stopifnot(is(basis,"BayesSpaceBasis"))
    if( !is.null(simulator) ) stopifnot(is.function(simulator))
    structure(list(density=density,
                   basis=basis,
                   paramSimulator=paramSimulator,
                   simulator=simulator,
                   hyperparameters=.hyperparameters),class="B2DistHyper")
}



#' Simulates a B2DistHyper distribution
#'
#'
#'
#' @param object an B2DistHyper object as returned from B2DistHyper
#' @param nsim the number of simulations
#' @param seed a seed eventually to be set before the simulation (use with care!)
#' @param density A function that can be called with a vectorized x of
#'     data values of the BayesSpace and parameters, to calculate a
#'     density.
#' @param basis A BayesSpaceBasis (or an object that can be converted
#'     to one) of BayesSpace containing the values of the distribution
#'     (of distributions).
#' @param paramSimulator A function to be called without arguments to
#'     simulate one set of parameters.
#' @param hyperparameters the list of parametrs to be given to the the
#'     parameterSimulator
#' @param simulator A function that can be called with an integer n
#'     and the parameters to create n random simulations.
#' @param ... futher unused parameters, required to fullfill the
#'     generic definition of simulate
#' @returns A [BDataDensities] following the given distribution.
#'
#' @seealso [stats::simulate]
#  @family B2DistHyper
#' @export
#' @examples
#' D = B2DistHyper(dbeta,B2Unif01Basis5,function() rexp(2),simulator=rbeta)
#' simulate(D,4,seed=47)
#'
simulate.B2DistHyper = function(object, nsim = 1, seed = NULL,
                                density=object$density,
                                basis=object$basis,
                                paramSimulator=object$paramSimulator,
                                simulator=object$simulator,
                                hyperparameters=object$hyperparameters,...
                                ) {
    if( !is.null(seed) ) {
    #    oState = .Random.seed
        set.seed(seed)
    }
    X = lapply(1:nsim,function(i) {
        theta=lapply(do.call(paramSimulator,hyperparameters),function(x) x)
        if( !is.null(simulator) ) {
            rand = function(n) do.call(simulator,c(list(n),theta))
        }
        f   = function(x) do.call(density,c(list(x),theta))
        f
    }
    )
    #if( !is.null(seed) ) {
      # .Random.seed <<- oState # it is not allowed by CRAN policies to change .Random.seed
    #}
    structure(BDataDensities(X,basis=object$basis),seed=seed)
}

#' A Gaussian Process distribution for probabilities
#'
#' @param mean something covertable to a single BDataLogDensity,
#'     providing the Bayes space expectation of the distribution, as
#'     well as BayesSpace and Support. Currently only finite continuous
#'     supports are possible. You might want to simulate on B2Unif01
#'     and transform for other Transformed spaces.
#' @param var a function of two arguments providing the covariance
#'     operator
#' @param nEval the number of equispaced discretization points used on
#'     the interval.
#' @param object the B2GaussianDistribution object
#' @param nsim the number of simulations to be computed
#' @param seed the optional seed, if used it will be reported as a
#'     seed attribution of the result (use with care)
#' @param ... unused arguments for generic use
#'
#' @returns B2GaussianDistribution: a B2GaussianDistribution object
#'     which can be simulated with, simulate: BDataLogDensities-Object
#'     with nsim cases with the same BayesSpace and BayesSpaceBasis as the mean.
#' @seealso [
#' @export
#'
#' @examples
#' m = BDataPDFs(list(dunif),basis=B2Unif01Basis7)
#' G1 = B2GaussianDistribution(m,function(x,y) exp(-(4*abs(x-y)^1.7)))
#' plot(simulate(G1,10,seed=47),col=2:11)
#'
B2GaussianDistribution = function(
                      mean=BDataLogDensities(list(function(x) rep(0,length(x))),basis=B2Unif01Basis7),
                      var=varKernel(\(x,y) exp(-abs(outer(x,y,"-"))),
                                    basis=getBayesSpaceBasis(mean)),
                      nEval=513) {
    stopifnot(numberOfCases(mean)==1)
    if( is(var,"varKernel") ) {
        # good we can use that
    } else if(is(var,"variationKernel")) {
        # good we can use that
    } else if(is(var,"function")) {
        var = varKernel(makeFunctionOuter(var),basis=getBayesSpaceBasis(mean))
    } else {
        stop("B2GaussianDistribution needs a varKernel, variationKernel, covariance function or 2-ary variogram as vari parameter")
    }
    structure(list(mean=mean,
                   var=var,
                   nEval=nEval),class="B2GaussianDistribution")
}

matSqrt = function(x) {
    result = with(svd(x),u%*%diag(sign(d)*sqrt(abs(d)))%*%t(v))
    SQR = result%*%t(result)
    D   = x-SQR
    if( any(!is.finite(D) | abs(D)>1E-6) ) {
        warning("matSqrt did not reproduce x")
        i = which(!is.finite(D) | abs(D)>1E-6)
        message("Mismatches at...")
        message(data.frame(i=i,col=col(x)[i],row=row(x)[i],x=x[i],rrt=D[i]))
    }
    result
}
#' Simulate a Gaussian Distribution on B2
#' @describeIn B2GaussianDistribution Simulaton method
#' @importFrom stats simulate
#' @importFrom stats rnorm
#' @seealso [stats::simulate]
#' @export
#' @method simulate B2GaussianDistribution
#'
#' @examples
#' GD = B2GaussianDistribution(mean=BDataPDFs(\(x)x,basis=B2Unif01Basis5),
#'                             var=\(x,y) exp(-abs(x-y)))
#' plot(simulate(GD,5))
#' GD = B2GaussianDistribution(mean=BDataPDFs(\(x) dnorm(x),basis=B2N01Basis5),
#'                             var=\(x,y) exp(-abs(x-y)))
#' plot(as.BDataPDFs(simulate(GD,5)))
#'
simulate.B2GaussianDistribution = function(object, nsim = 1, seed = NULL, ..., nEval=object$nEval) {
    if( !is.null(seed) ) {
    #    oState = .Random.seed
        set.seed(seed)
    }
    mean = as.BDataLogDensities(object$mean)
    basis=getBayesSpaceBasis(mean)
    D=getDiscretizer(basis)
    ym = mean[[1]](D$all)
    C = object$var(D$all,D$all)
    reportIfNot(ncol(C)==nrow(C),dimC=dim(C))
    reportIfNot(ncol(C)==length(D$all),dimC=dim(C),lenDall=length(D$all))
    reportIfNot(all(is.finite(C)),nf=which(!is.finite(C)))
    reportIfNot(all(abs(C)<1.1),nf=which(!(abs(C)<1.1)))
    SMat = matSqrt(C)
    fs   = replicate(nsim,
                     DiscretizedFunction(D,c(ym+SMat%*%stats::rnorm(ncol(SMat))))
                     )
    #if( !is.null(seed) ) {
      # .Random.seed <<- oState # it is not allowed by CRAN policies to change .Random.seed
    #}
    BDataLogDensities(fs,basis=basis,seed=seed)
}