There are a class of functions called ‘probability generating functions’ and ‘moment generating functions’ which are useful for describing various qualities about both discrete and continuous functions.  These functions are sometimes also referred to as the ‘expectation’ of the function.   The simplest of these is simply the first moment, E(X) and can be thought of as the mean value of the function over it’s defined range.   It is given by:  E(X) = ? x*f(x) dx.    The second moment is: E(X) = ? x²*f(x) dx.  The discrete version of these can be found by simply summing over the defined range.   Moment generating functions are generally defined as m_x(t) = E(e^sX) where s ? R¹.

# Calculate the moment of a random variable.  Currently only functional
# for a single variable.
# eg. probDensFun <- function(x, n=20, theta=0.5) {  choose(n,x)*((1-theta)^(n-x))*theta^x  }  # binomial
#     probDensFun <- function(x) { y=x %in% 1:6; y=y*1/6 }  # dice roll
#     probDensFun <- function(x) { x^2 }                    # x squared
#     probGenFun  <- function(x) {  x  }                # get basic expectation (ie mean)
#     probGenFun  <- function(x, s=0) {  exp(s*x)  }    # moment generating function
#     xSeq        <- 1:20
#     expectation(probDensFun, probGenFun, xSeq)
expectation <- function(probDensFun=function(x){},
                        probGenFun=function(x){ x },   # default = expectation
                        xSeq=seq(0,1, length.out=10),
                        discrete=FALSE) { 

  require(MASS)
  # If the density function is continuous, then use sophisticated numerical
  # integration techniques to estimate the area under the curve.
  if (!discrete) {
    require(cubature)
    f        <- function(x) probDensFun(x) * probGenFun(x)
    result   <- adaptIntegrate(f, lowerLimit=min(xSeq), upperLimit=max(xSeq))
    result$f <- paste(c(deparse(probGenFun, control="all"), "*",
                        deparse(probDensFun, control="all")),
                      collapse=" ")
    result$x         <- range(xSeq)
    result$fraction  <- fractions(result$integral)
    return(result)

  # Otherwise, this is a discrete problem so the program will multiply the
  # density function times the probability generating function at the
  # specified points along x.
  } else {
    result          <- list()
    dx              <- diff(xSeq);  dx <- c(dx, dx[length(dx)])
    result$integral <- sum(probDensFun(xSeq) * probGenFun(xSeq))
    result$f        <- paste(c(deparse(probGenFun, control="all"), "*",
                        deparse(probDensFun, control="all")),
                      collapse=" ")
    result$x        <- xSeq
    result$fraction <- fractions(result$integral)
    return(result)
  }
}