R/integration.R

In DistributionalDDA: Distributional Data Analysis in Bayes Spaces

## Creates a equispace uniform grid sample with half grid distance to ends
makeUnifGrid = function(n=2^16,a=0,b=1) {
    eps = (b-a)/n
    seq(a+eps/2,b-eps/2,length.out=n)
}


## Closes a weights vector
.closeWeights = function(w) w/sum(w)

## Computes mean from closed weights
.wMean = function(X,w) {
    if( is.matrix(X) ) {
        c(w %*% X)
    } else {
        sum(w*X)
    }
}

## Centers a dataset from closed weights
.wCenter = function(X,w) {
    if( is.matrix(X) ) {
        t(t(X)-c(w %*% X))
    } else {
        X-sum(w*X)
    }
}

.centerAround = function(X,m) {
    if(is.matrix(X)) {
        X-structure(rep(m,each=nrow(X)),dim=dim(X))
    } else {
        X-m
    }
}

## Computes the covariance of a dataset from closed weights
.wCov = function(X,Y,w) {
    X = .wCenter(X,w)
    Y = .wCenter(Y,w)
    .wE2(X,Y,w)
}

## Computes the uncentered covariance of a dataset from closed weights
.wE2 = function(X,Y,w) {
    if( is.matrix(X) ) {
        if( is.matrix(Y) ) {
            (t(X*w)%*%Y)
        } else {
            (t(X*w)%*%Y)
        }
    } else {
        if( is.matrix(Y) ) {
            t(X*c(w)) %*% Y
        } else {
            c((X*c(w)) %*% Y)
        }
    }
}

.wdiagE2 = function(X,w) {
    if(is.matrix(X)) {
        colSums(w*(X^2))
    } else {
        sum(w*X^2)
    }
}

## Computes Expectation using a sample
### @param X the sample (weighed or unweighted)
### @param f the real funcition on the dataset X to compute the values
### @param further arguments to f
### @return a numeric with the (column wise expectation)

.expectationUsingSample <- function(X,f,...) .averageForSample(X,f(X,...))

.expectationUsingChunk <- function(measure,f,...,junksize) {
    x <- simulate(measure,nsim=junksize)
    X <- f(x,...)
    .meanstats1ForSample(x,X)
}

.covUsingChunk <- function(measure,X,Y,...,junksize) {
    x <- simulate(measure,nsim=junksize)
    if( is.null(Y) )
        .varStats1ForSample(x,X(x,...))
    else
        .covStats1ForSample(x,X(x,...),Y(x,...))
}


.expectationUsingSim <- function(measure,f,...,tol=1E-3,
                                 junksize=1000000,
                                 nmax=1000000) {
    N = 0
    i=1
    A <- .expectationUsingChunk(measure,f,...,junksize=junksize)
    nsim=junksize
    m = A$m
    mv = A$mv
    N = N+A$N
    if(tol<=1E-16) tol=1E-16
    for(k in 1:nmax) {
        effTol = pmax(tol,abs(m)*tol)

        if( all(sqrt(mv) * 2 <= effTol) ) {
            if( i>1 )
                message(i,". MonteCarlo step: samples=",N," done=",min((effTol^2)/(4*mv)),"\nDone\n")
            return(structure(m,sd=sqrt(mv),nsim=N))
        } else {
            if(i==1)
              message("Iterative Simulation started, for tol=",tol,"\n")
            message(i,". MonteCarlo step: samples=",N," done=",min((effTol^2)/(4*mv)),"\n")
        }

        A = .expectationUsingChunk(measure,f,...,junksize=junksize)
        i=i+1
        Nnew = N+A$N
        m = m *(N/Nnew)+ A$m*(A$N/Nnew)
        mv = (mv*((N/Nnew)^2)+A$mv*((A$N/Nnew)^2))
        N = Nnew
    }
    warning("Expectation computation could not reach expected precisions in ",
            nmax, " cycles")
    return(structure(m,sd=sqrt(mv),nsim=N))
}
## Computes a mean and covariance.
## The tolerance unfortunately is only for the mean.
## If Y is not given, it defaults to X but with speedups

.covUsingSim <- function(measure,X,Y=NULL,...,tol=1E-3,
                                 junksize=1000000,
                         nmax=1000000,
                         weightFunction=getWeights) {
    N = 0
    i=1
    A <- .covUsingChunk(measure,X,Y,...,junksize=junksize)
    nsim=junksize
    mX = A$mX
    mY = A$mY
    mvX = A$mvX
    mvY = A$mvY
    C  = A$C
    N = N+A$N
    if(tol<=1E-16) tol=1E-16
    while(TRUE) {
        effTolX = pmax(tol,abs(mX)*tol)
        effTolY = pmax(tol,abs(mY)*tol)
        if( all(sqrt(mvX) * 2 <= effTolX) && all(sqrt(mvY) * 2 <= effTolY) ) {
            if( i>1 )
              message(i,". MonteCarlo step: samples=",N," done=",min((effTolX^2)/(4*mvX),(effTolY^2)/(4*mvY)),"\nDone\n")
            return(structure(mX=mX,mY=mY,sdX=sqrt(mvX),sdY=sqrt(mvY),C=C,nsim=N))
        } else {
            if(i==1)
              message("Iterative Simulation started, for tol=",tol,"\n")
            message(i,". MonteCarlo step: samples=",N," done=",min((effTolX^2)/(4*mvX),(effTolY^2)/(4*mvY)),"\nDone\n")
        }

        A = .covUsingChunk(measure,X,Y,...,junksize=junksize)
        i=i+1
        Nnew = N+A$N
        m = m *(N/Nnew)+ A$m*(A$N/Nnew)
        C = C *(N/Nnew)+ A$C*(A$N/Nnew)
        mvX = (mvX*((N/Nnew)^2)+A$mvX*((A$N/Nnew)^2))
        mvY = (mvY*((N/Nnew)^2)+A$mvY*((A$N/Nnew)^2))
        N = Nnew
    }
}


## Averages over a sample generated from X (which might have weights)
.averageForSample <- function(X,Y,weightFunction=getWeights) {
    N    <- numberOfCases(X)
    hasW <-   hasWeights(X)
    if( length(Y) == N ) {
        if(hasW)
            stats::weighted.mean(Y,c(weightFunction(X)))
        else
            mean(Y)
    } else {
        reportIfNot(nrow(Y)==N,dim(Y),N,"Expected data driven structure")
        if(hasW) {
            w = c(weightFunction(X))
            reportIfNot(length(w)==N,length(w),N,"Expected correct number of weights")
            c((w/sum(w)) %*% Y)
        } else
            colMeans(Y)
    }
}
## Computes weighted mean and invidivual variance for uni- or multivariate data
## Returns
##  m  the vector of means
##  mv the variance of the means
.meanstats1ForSample <- function(X,Y,weightFunction=getWeights,close=FALSE) {
    N    <- numberOfCases(X)
    reportIfNot(nCases(Y)==N,dim(X), length(X),dim(Y), length(Y),"X does not fit to Y")
    w = c(weightFunction(X))
    if( close ) {
        w=w/sum(w)
    }
    if( length(Y) == N ) {
        Z = w*Y
        m = sum(Z)
        list(m=m,mv=var(Z)*N,N=N)
    } else {
        Z  = w*Y
        m  = colSums(Z)
        mv = colSums((Z-outer(w,m))^2)*(N/(N-1))
        list(m=m,mv=mv,N=N)
    }
}

.covStats1ForSample = function(x,X,Y,weightFunction=getWeights,close=FALSE) {
    N    <- numberOfCases(x)
    reportIfNot(numberOfCases(X)==N || numberOfCases(Y)==N,
                dim(x),dim(X), dim(Y), length(x), length(X),length(Y),"x does not fit to X or Y")
    w = c(weightFunction(x))
    sw= sum(w)
    if( close ) {
        w=w/sw
        sw=1
    }
    if( length(X)==N && length(Y)==N ) {
        mX = sum(w*X)
        mY = sum(w*Y)
        vX = var(w*X)
        vY = var(w*Y)
        C = sum(w*(X-mX/sw)*(Y-mY/sw))
        list(mX=mX,mY=mY,mvX=vX*N,mvY=vY*N,C=C,N=N)
    } else {
        X = as.matrix(X)
        Y = as.matrix(Y)
        mX = colSums(w*X)
        mY = colSums(w*Y)
        MX=NCOL(X)
        MY=NCOL(Y)
        vX = apply(X,2,function(x) var(w*x))
        vY = apply(Y,2,function(y) var(w*y))
        C = .MatYApply(1:MY,function(j) .MatXApply(1:MX,function(i)  sum(w*(X[,i]-mX[i]/sw)*(Y[,j]-mY[j]/sw))))
        list(mX=mX,mY=mY,mvX=vX*N,mvY=vY*N,C=C,N=N)
   }
}

.varStats1ForSample = function(x,X,weightFunction=getWeights,close=FALSE) {
    N    <- numberOfCases(x)
    reportIfNot(numberOfCases(X)==N ,
                dim(x),dim(X),length(x),length(X),"x does not fit to X")
    w = c(weightFunction(x))
    sw= sum(w)
    if( close ) {
        w=w/sw
        sw=1
    }
    if( length(X)==N ) {
        mX = sum(w*X)
        vX = var(w*X)
        C = sum(w*(X-mX/sw)^2)
        list(mX=mX,mY=mX,mvX=vX*N,mvY=vX*N,C=C,N=N)
    } else {
        X = as.matrix(X)
        mX = colSums(w*X)
        MX=NCOL(X)
        vX = apply(X,2,function(x) var(w*x))
        C = .MatYApply(1:MX,function(j) .MatXApply(1:MX,function(i)  sum(w*(X[,i]-mX[i]/sw)*(X[,j]-mX[j]/sw))))
        list(mX=mX,mX=mX,mvX=vX*N,mvY=vX*N,C=C,N=N)
   }
}


#' Computes a covariance with respect to a measure. The function can
#' be considered internal implementation details of the package, but
#' needs to be exported due to its generic nature.
#'
#' The implementation depends on the measure
#'
#' @param measure typically a ReferenceMeasure
#' @param X a random variable (or a list of such) on the domain of the
#'     measure represented as an R function. It returns multiple
#'     values as columns of a matrix.
#' @param Y a random variable (or a list of) on the domain of the
#'     measure represented as an R function.
#' @param ... possibly further implementation dependent arguments for
#'     the computational method
#' @param close if false the underlying measure is not assumed to be
#'     closed, and the function makes it possible to computed unclosed
#'     covariance integrals for unclosed reference measures. This is
#'     the default.
#'
#' @returns the scalar products (covariances) of `X` and `Y` with
#'     respect to the reference measure. For the default
#'     ReferenceMeasure method these are obtained from [stats::cov()]
#'     applied to `X` and `Y` evaluated on the measure's sample. The
#'     result is a numeric scalar, vector, or matrix depending on the
#'     number of columns returned by `X` and `Y` (one entry per pair of
#'     output components); when `X` or `Y` is a list, the results are
#'     combined columnwise.
#'
#' @export
#' @importFrom stats weighted.mean
#'
#' @examples
#' CovarianceOperator(StdReferenceMeasureUnif01,
#'                    function(x) cbind(x,-x),
#'                    function(x) cbind(x,-x,2*x+1)
#' )
#'
CovarianceOperator = function(measure,X,Y,...) UseMethod("CovarianceOperator",measure)
#' @export
#' @describeIn CovarianceOperator method for ReferenceMeasures using a
#'     MonteCarlo integration based on the sample provided by
#'     getSample(sample)
#' @param method a method to be selected. Typical methods are
#'     "sample", "MonteCarlo" and "integrate", but only sample is
#'     availabel for all reference measures. "sample" uses the
#'     prestored sample, "MonteCarlo" uses the simulated function of
#'     the reference measure, and "integrate" uses numerical
#'     integration. Typically the prestored integration sample
#'     performs well with derivatives measures derived
#'     ReferenceMeasureUniform, ReferenceMeasusureTransformed, where
#'     it proves a dense grid in the area of measure, and
#'     ReferenceMeasure discrete where it uses precise computation
#'     using a weighted sample. For Cartesian products involving more
#'     than one continues measure due to the curse of dimensionality
#'     typically only MonteCarlo is expected to work reasonably well.
#' @param tol the absolute or relative tolerance (whichever is larger).
#'     accepted for methods with tolerance control.

#' @importFrom stats cov
#' @method CovarianceOperator ReferenceMeasure
CovarianceOperator.ReferenceMeasure = function(measure,X,Y,...,method=c("default","sample","MonteCarlo","integration"),tol=1E-2,close=FALSE) {
    method <- selectIntegrationMethod(measure,match.arg(method))
    if( is.list(Y) ) {
        .MatYApply(Y,function(Y1) CovarianceOperator(measure,X,as.function(Y1),method=method,tol=tol,close=close))
    } else if( is.list(X) ) {
        .MatXApply(X,function(X1) CovarianceOperator(measure,as.function(X1),Y,method=method,tol=tol,close=close))
    } else  {
        if( method=="sample" ) {
            x = getSample(measure)
            .covStats1ForSample(x,X(x),Y(x),close=close)$C
        } else {
            if( method!="MonteCarlo" ) {
                warning("Unimplemented method  ",method," in CovarianceOperator : using MonteCarlo")
            }
            .covUsingSim(measure,X,Y,...,tol=tol,close=close)$C
        }
    }
}

.MatXApply = function(X,FUN,...) {
    ## Cov(X,Y)          X FUN
    ## FUN Scalar        N 1  => N x 1
    ## FUN Vector        N M  => N x M
    ## FUN Spaltenmatrix N M  => sollte nicht passieren
    ## FUN ZeilenMatrix  N M  => N x M
    ## FUN Scalar        1 1  => 1 x 1
    ## FUN Vector        1 M  => 1 x M
    ## FUN Spaltenmatrix 1 M  => sollte nicht passieren
    ## FUN ZeilenMatrix  1 M  => 1 x M
    C = t(matrix(sapply(X,FUN,...),ncol=length(X)))
}

.MatYApply = function(Y,FUN,...) {
    ## Cov(X,Y)        FUN Y
    ## FUN Vector        N 1  => N x 1
    ## FUN Vector        N M  => N x M
    ## FUN Spaltenmatrix N M  => N x M
    ## FUN ZeilenMatrix  N M  => sollte nicht passieren
    ## FUN Scalar        1 1  => 1 x 1
    ## FUN Scalar        1 M  => 1 x M
    ## FUN Scalar        1 M  => 1 x M
    ## FUN Scalar        1 M  => 1 x M
    C = matrix(sapply(Y,FUN,...),ncol=length(Y))
}



#' Computes mean and variance for weighted reference measures
#'
#' The function computes the expectation and variance of 'X' for the
#' probability distribution given by the density f with respect to the
#' measure.  The method supports lists even if the individual measure
#' does not.
#'
#' The default method for reference measures used [getSample] and
#' Monte Carlo integration with respect to this sample to perform the
#' computation. More advanced methods e.g. using numerical integration
#' might be implemented for certain subclasses.
#'
#' @param measure should be a reference measure. Other classes can
#'     implement the method.
#' @param X a function (or list of functions) providing the statistics
#' @param f a function providing the density with respect to X
#' @param method a method to be selected. Typical methods are
#'     "sample", "MonteCarlo" and "integrate", but only sample is
#'     availabel for all reference measures. "sample" uses the
#'     prestored sample, "MonteCarlo" uses the simulated function of
#'     the reference measure, and "integrate" uses numerical
#'     integration. Typically the prestored integration sample
#'     performs well with derivatives measures derived
#'     ReferenceMeasureUniform, ReferenceMeasusureTransformed, where
#'     it proves a dense grid in the area of measure, and
#'     ReferenceMeasure discrete where it uses precise computation
#'     using a weighted sample. For Cartesian products involving more
#'     than one continues measure due to the curse of dimensionality
#'     typically only MonteCarlo is expected to work reasonably well.
#' @param tol the absolute or relative tolerance (whichever is larger).
#'     accepted for methods with tolerance control.
#'
#' @returns a list with the elements E for the computed expectation
#'     and V for the computed variance.
#' @seealso [getSample]
#' @export
#' @examples
#' weightedStats(ReferenceMeasureUniform(0,1),function(x) x, function(x) x^2)
#' weightedStats(ReferenceMeasureUniform(0,1),as.function(B2Unif01Basis5), function(x) x^2)
weightedStats   = function(measure,X,f,method=c("default","sample","MonteCarlo","integration"),tol=1E-2) {
    method <- selectIntegrationMethod(measure,match.arg(method))
    if( is.list(X) )
        as.data.frame(t(sapply(X,function(x) weightedStats(measure,X=as.function(x),f=as.function(f),method=method,tol=tol))))
    else
        UseMethod("weightedStats",measure)
}

#' @describeIn weightedStats method for ReferenceMeasure
#' @param ... further arguments to integration methods e.g. tol.
#' @export
#'
weightedStats.ReferenceMeasure = function(measure,X,f,method=c("default","sample","MonteCarlo","integration"),...) {
    method <- selectIntegrationMethod(measure,match.arg(method))
    if( method=="sample" ) {
        # x = getSample(measure)
        # w = getWeights(x)
        # XX = as.matrix(X(x))
        # wt=c(f(x))*w
        # wt=wt/sum(wt)
        # reportIfNot(nrow(XX)==length(wt),dimXX=dim(XX),lengthwt=length(wt),dimwt=dim(wt))
        # val = stats::cov.wt(XX,wt=wt,method="ML")
        # list(E=val$center,V=val$cov)
        x = getSample(measure)
        with(.varStats1ForSample(x,X(x),weightFunction=\(x) getWeights(x)*f(x),
                                 close=TRUE),
             list(E=mX,V=C))
    } else {
        if(method!="MonteCarlo") warning("unsupported integration method:",method, "using MonteCarlo.")
        with(.covUsingSim(measure,X,Y=NULL,...,
                          weightFunction=function(x) getWeights(x)*f(x),
                          close=TRUE),
             list(E=mX,V=C))
    }
}