# R/symmoments.R In symmoments: Symbolic central and noncentral moments of the multivariate normal distribution

#### Defines functions toSquare

multmoments <-
function (moment, current.matrix, current.cell, moment.rep, row_col)
{
# A recursive function to compute the representation of a multivariate moment
#    using upper - triangular matrices
#
# moment: a vector of integers representing the moment
#     eg, c(3, 2, 4) for a 3-dimensional normal vector
#     corresponding to the moment E[(X1^3)(X2^2)(X3^4)]

# current.matrix: input/output matrix under consideration in recursion
#     this is an upper-triangular integer matrix (l(i, j))
#     organized by row, with   lmom * (lmom + 1) / 2   total elements
#
# current.cell: input/output cell in current matrix under consideration in recursion

# moment.rep: input/output current set of representations
#     function adds each satisfying matrix to moment.rep

# row_col: 2x(lmom * (lmom + 1) / 2 matrix giving rows and columns for square matrix
#     for each cell so that they don't have to be calculated each time

#  algorithm:
#  loop through cell values from 0 to min(moments[row] - rowsum, moments[col] - colsum)
#    if the this new matrix satisfies moment criterion,
#       add to moment.rep and return
#    if the current matrix is too great in any dimension,
#       return
#    if the current matrix is < moment in any dimension, and
#       at most moment(i) for all other indexs i, continue

summomentmatrix <- function (moment.matrix)
{
#  compute the row/col sums of moment.matrix
#  uses Matrix package
#  construct a square upper diagonal matrix from moment.matrix

#  length of moment based on representation:
lmom <- (sqrt(8 * length(moment.matrix) + 1) - 1) / 2
#  make a square matrix for use with row and column sum functions
tempmatrix <- matrix(c(rep(0, lmom^2)), nrow=lmom)
tempmatrix[1, ] <- moment.matrix[1:lmom]
endrow <- lmom
if (lmom > 1)
{
for (irow in 2:lmom)
{startrow <- endrow + 1
endrow <- endrow + (lmom - irow + 1)
tempmatrix[irow, ] <- cbind(c(rep(0, (irow - 1)), moment.matrix[startrow:endrow])) }

summomentmatrix <- colSums(tempmatrix) + rowSums(tempmatrix)
}
if (lmom == 1){summomentmatrix <- 2 * sum(moment.matrix)}

return(summomentmatrix)}

lmom <- length(moment)
totcells <- lmom * (lmom + 1) / 2   # total cells in a moment representation

thisrow <- row_col[1, current.cell]
thiscol <- row_col[2, current.cell]

moment.row <- moment[thisrow]
moment.col <- moment[thiscol]
rowcells <- row_col[1, row_col[1, ] == thisrow]
colcells <- row_col[2, row_col[2, ] == thiscol]

#  determine sums for row and columm
#  maxvalue is the largest value that can be added to a
#    a row/column sum and still be no more than criterion
#  maxvalue will be the minimum of the moments minus these sums

rowsum <- summomentmatrix(current.matrix)[thisrow]
colsum <- summomentmatrix(current.matrix)[thiscol]

maxvalue <- min(moment[thisrow] - rowsum, moment[thiscol] - colsum)

for (ivalue in 0:maxvalue)
{         # for:loop through all possible values for cell

current.matrix[current.cell] <- ivalue
current.sum <- summomentmatrix(current.matrix)
# determine current sum for this matrix

#  if matrix fulfills criterion, add to reps and return
if (sum(moment == current.sum) == lmom)
{moment.rep <- rbind(moment.rep, current.matrix)
return(moment.rep)}

#   if the sum is too large, return because any other sum will also be too large
if (sum(current.sum > moment) > 0){return(moment.rep)}

#   if at least one term in current sum is smaller than moment,
#       go down one cell unless this is the last cell

if ((sum(current.sum == moment) < lmom) & (current.cell != totcells))
{cc <- current.cell + 1
# recursive step:
moment.rep <- multmoments(moment, current.matrix, cc, moment.rep, row_col)}

}        # end of for

#   note: the three conditions above are exhaustive except for last cell
#      In that case there is nothing to return
return(moment.rep)
}

callmultmoments <-
function (moment)
{
#  function to compute the representation of a multivariate moment
#
#  moment: a vector of integers representing the moment
#          eg, c(3, 2, 4) for a 3-dimensional normal vector
#          corresponding to the moment (E[X1^3)(X2^2)(X3^4)]
#  sum of exponents must be even; otherwise moment is 0

#  returns a list of 3 components:
#    1: $moment - the input moment vector # 2:$representation  -  a matrix containing the representation
#       in terms of upper-triangular matrices
#    3: $coefficients - the coefficients corresponding to the rows of the representation # if sum odd, returns -1 and prints "Sum of powers is odd. Moment is 0." # if any component is negative, returns -2 and prints "All components of the moment must be non-negative." # if any component is not an integer, returns -3 and prints "All components of the moment must be integers." subblank <- function (inputstring) { # remove blanks from a character string outputstring <- sub(" ", "", inputstring) if (outputstring == inputstring) { return(outputstring) } if (outputstring != inputstring) { outputstring <- subblank(outputstring) } return(outputstring) } nestedreps <- function (input.vector, inner.rep, outer.rep) { # replicates input.vector, first by inner.rep, then by outer.rep # input.vector: vector to replicate # inner.rep: count of replicates of elements # outer.rep: count of replicates of resulting vector temp.vector <- NULL for (ivec in 1:length(input.vector)) {temp.vector <- c(temp.vector, rep(input.vector[ivec], inner.rep))} total.vector <- rep(temp.vector, outer.rep) return(total.vector)} mrepnames <- function (ndim) { # get colnames for a representation combs <- expand.grid(1:ndim, 1:ndim) char.combs <- paste("S(", combs[, 2], ", ", combs[, 1], ")") for (elem in 1:(ndim^2)) {char.combs[elem] <- subblank(char.combs[elem])} m <- sort(matrix((1:ndim^2), nrow=ndim, byrow=TRUE)[ !lower.tri(matrix((1:ndim^2), nrow=ndim, byrow=TRUE))]) return(char.combs[m]) } lmom <- length(moment) lrep <- (lmom^2 + lmom) / 2 # length of representation using upper-triangular matrices moment.rep <- matrix(rep(0, lrep), nrow=1) # initial null representation to be augmented if (sum(trunc(moment) == moment) < length(moment)) {print("All components of the moment must be integers.") return( -3)} if (sum(moment<0) > 0) {print("All components of the moment must be non-negative.") return( -2)} if (trunc(sum(moment) / 2) != sum(moment) / 2) {print("Sum of powers is odd. Moment is 0.") return( -1)} icells <- matrix(rep(0, lmom^2), nrow=lmom) m <- matrix(rep(1, lmom^2), nrow=lmom) limits <- c((lmom:1)%*%(m * !(lower.tri(m)))) # sum of row lengths: lmom, lmom + lmom - 1, ... , for use in row_col row_col <- matrix(rep(0, (2 * lrep)), nrow=2) # 2x(nm * (nm + 1) / 2 matrix giving rows and columns for each cell # so that they don't have to be calculated each time in multmoment for (icell in (1:lrep)) { row_col[1, icell] <- min((1:lmom)[icell<=limits]) if (row_col[1, icell] == 1){row_col[2, icell] <- icell} if (row_col[1, icell]>1){row_col[2, icell] <- icell - limits[row_col[1, icell] - 1] + row_col[1, icell] - 1 } } # initial current.matrix and current.cell current.matrix <- c(rep(0, lrep)) current.cell <- 1 # call recursive function to determine upper-triangular representations moment.rep <- multmoments(moment, current.matrix, current.cell, moment.rep, row_col) if (dim(moment.rep)[1] == 2) # get rid of initial 0 matrix representation { moment.rep <- matrix(moment.rep[2, ], nrow=1) } if (dim(moment.rep)[1] > 2) {moment.rep <- moment.rep[2:dim(moment.rep)[1], ]} rownames(moment.rep) <- 1:(dim(moment.rep)[1]) ################################################################## # now determine coefficients for upper-triangular representations l.representation <- moment.rep lmom <- length(moment) nrep = dim(l.representation)[1] totlength <- lmom^2 rep.coefficients <- c(1:nrep) # coefficients corresponding to nrep representations # multiplier for all terms overallcoeff <- ((1 / 2)^(sum(moment) / 2)) * prod(factorial(moment)) / factorial((sum(moment) / 2)) for (irep in 1:(dim(l.representation)[1])) { # loop through all matrices thisrep <- l.representation[irep, ] # determine the coefficient for each term based on switching equivalent terms # "base" gives the number of switches that can be made to each element of the l-matrix # diagonal elements are not switchable, but are included to allow subtraction below base <- c(rep(1, lmom * (lmom + 1) / 2)) base[1] <- 1 # first diagonal element - not switchable totreps <- 1 # total number of transpostions # if there is only one element, it must be the diagonal, so is not switchable - skip if (lmom > 1){ base[1] <- 1 # first diagonal element for (cell in 2:length(base)) { icol = row_col[1, cell] # determine if diagonal element irow = row_col[2, cell] if (irow == icol){base[cell] <- 1} # diagonal - not switchable if (icol != irow) {totreps <- totreps * (1 + thisrep[cell]) base[cell] <- 1 + thisrep[cell]} } # done with computing base and total transpositions (totreps) } mcoeff <- 1 # sum of multinomial coefficients if (totreps > 1){ # baserep will represent the lower diagonal (including diagonal) # of the augmented matrices baserep <- matrix(rep(0, totreps * length(base)), nrow=totreps) basegt1 <- base[base>1] nbase <- 0 for (ibase in 1:length(base)) {if (base[ibase] > 1) {nbase = nbase + 1 if (nbase == 1){baserep[, ibase] <- nestedreps(c(0:(basegt1[nbase] - 1)), 1, totreps / prod(basegt1[1:nbase])) } if (nbase > 1) {baserep[, ibase] <- nestedreps(c(0:(basegt1[nbase] - 1)), prod(basegt1[1:(nbase - 1)]), totreps / prod(basegt1[1:nbase])) } } } # now go through each transposition if ( !is.na(totreps) & totreps != 1) { mcoeff <- 0 for (jrep in (1:totreps)) # check each transposition {newrep <- baserep[jrep, ] # added lower diagonal elements addrep <- sort(newrep, decreasing=TRUE)[1:(lmom * (lmom - 1) / 2)] fulnrep <- c((thisrep - newrep), addrep) thiscoeff <- ((length(fulnrep))^sum(fulnrep)) * dmultinom(x=fulnrep, prob=rep(1.0, length(fulnrep))) mcoeff <- mcoeff + thiscoeff # the multinomial coefficient is obtained from the multinomial distribution # multiply by an appropriate power to get rid of probability } } } if (is.na(totreps)){mcoeff <- 1} if (totreps == 1) {mcoeff <- (length(thisrep))^sum(thisrep) * dmultinom(x=thisrep, prob=rep(1.0, length(thisrep)))} # determine full coefficient - round because all coefficients should be integers # (Note - this statement has not been proved) rep.coefficients[irep] <- round(overallcoeff * mcoeff) cell <- 0 for (irow in (1:length(moment))) { for (icol in (irow:length(moment))) {cell <- cell + 1 # exponent of term expo <- l.representation[irep, cell] } } } output <- list(moment, moment.rep, rep.coefficients) names(output) <- c('moment', 'representation', 'coefficients') colnames(output$representation) <- mrepnames(length(moment))
names(output$coefficients) <- paste("rep", (1:length(output$coefficients)))
class(output) <- "moment"
return(output)}

toLatex.moment <-
function (object, ...)
{
#  build latex code for the l-matrix representation of the moment
#  object is the representation of the l-matrices for moment
#  each row is such an l-matrix

# object: list from callmultmoment (class moment)
#     with first component the moment itself,
#     the second component the set of upper-triangular
#          representations of the moment,
#     and third component, their correpsonding coefficients

#  note that Latex backslashes are doubled to allow writing to file
#      with writeLines

subblank <- function (inputstring)
{
# remove blanks from a character string
outputstring <- sub(" ", "", inputstring)
if (outputstring == inputstring) {
return(outputstring)
}
if (outputstring != inputstring) {
outputstring <- subblank(outputstring)
}
return(outputstring)
}
sortmatrix <- function (xmat)
{
# sort matrix xmat by successive columns
# output is the index vector, not the sorted matrix

ncol <- dim(xmat)[2]

sortstring <- "order("
for (icol in 1:ncol)
{
sortstring <- paste(sortstring, "xmat[, ", icol, "]")
if (icol != ncol){sortstring <- paste(sortstring, ", ")}
}
sortstring <- paste(sortstring, ", decreasing=TRUE)")
matindex <- eval(parse(text=sortstring))

return(matindex)}

moment.fullrep <- object

maxchars.latex <- 300
# maximum latex characters to print on one line
# this includes formatting characters (subscript notation, etc)

#  extract the components for convenience
moment <- moment.fullrep[[1]]
moment.rep <- moment.fullrep[[2]][sortmatrix(moment.fullrep[[2]]), ]
#       latex representation will be in sorted order
coefficients.rep <- moment.fullrep[[3]][sortmatrix(moment.fullrep[[2]])]

doubquote <- subblank(paste("\\", "\\"))

if ( !is.matrix(moment.rep)){moment.rep <- matrix(moment.rep, nrow=1)}
numrep <- dim(moment.rep)[1]

latex.moment <- rep(" ", numrep + 1)

#   write the left hand side, ie, E[X1 ... Xn] =

latex.moment[1] <- "E["

for (imoment in(1:(length(moment))) )
{latex.moment[1] <- paste(latex.moment[1],
"X_{", imoment, "}^{",  moment[imoment],  "}")  }

latex.moment[1] <- paste(latex.moment[1], "] =", doubquote)
latex.moment[1] <- subblank(latex.moment[1])

#  write the right hand side, ie, the set of terms

if (sum(moment==0) == length(moment))
{latex.moment[2] <- "1"
return(latex.moment)}

totchars <- 0        # used with totchars.latex
for (irep in (1:numrep))
{
mcoeff <- coefficients.rep[irep]
thisrep <- moment.rep[irep, ]

if (mcoeff != 1){latex.moment[irep + 1] <- as.character(mcoeff)}
#              omit coefficient if it is 1
cell <- 0
for (irow in (1:length(moment)))
{
for (icol in (irow:length(moment)))
{cell <- cell + 1
#   exponent of term, that is, the "l" value
exponent <- moment.rep[irep, cell]

#   if exponent is 1, omit it as obvious
if (exponent == 1){latex.moment[irep + 1] <- paste(latex.moment[irep + 1],
"\\sigma_{ ",  irow,  ", ",  icol,  "}")}

if (exponent > 1){latex.moment[irep + 1] <- paste(latex.moment[irep + 1],
"\\sigma_{ ",  irow,  ", ",  icol, "}^{",  exponent, "} ")}
}  # end of for
}  # end of for
totchars <- totchars + nchar(latex.moment[irep + 1])
if (irep < numrep){latex.moment[irep + 1] <- paste(latex.moment[irep + 1], " + ")}
if (totchars > maxchars.latex)
{totchars <- 0
latex.moment[irep + 1] <- paste(latex.moment[irep + 1], doubquote)}
latex.moment[irep + 1] <-  subblank(latex.moment[irep + 1])
}  # end of for

return(latex.moment)}

simulate.moment <-
function(object, nsim, seed=NULL, Mean, Sigma, ...){

# function: method to calculate moment of the multivariate normal distribution
#           using Monte-Carlo integration (Rizzo, 2008)
# object is an object of class moment
# nsim is the number of samples to generate
# seed is the seed for the random number generator
# Mean is the mean of the (X1, ..., Xn)
# Sigma is the variance-covariance of (X1^k1, ..., Xn^kn), dimension nXn

# requires package mvtnorm for function rmvnorm

moment.fullrep <- object
if (is.numeric(seed)){set.seed(seed)}

if (class(moment.fullrep) == "moment"){thismoment <- moment.fullrep$moment} if (class(moment.fullrep) != "moment") {print("moment must be of class 'moment'") return(-1)} ndim <- length(thismoment) sample <- rmvnorm(n=nsim, mean=Mean, sigma=matrix(Sigma,nrow=length(Mean))) exponents <- matrix(rep(thismoment, nsim), nrow=nsim, byrow=TRUE) powers <- sample^exponents prods <- rep(1, nsim) # calculate product of powers of Xs for (icol in (1:ndim)) {prods <- prods * powers[, icol]} moment.value <- mean(prods) return(moment.value)} evaluate <- function(object, sigma) {UseMethod("evaluate",object)} evaluate.moment <- function (object, sigma) { # evaluate the moment using the representation from callmultmoment # at the upper-triangular value of sigma # object from callmultmoment # list with first component the moment itself # the second component the set of upper-triangular # matrices representing the moment # and third component, their corresponding coefficients # # sigma is the upper-triangular matrix of covariance terms # at which the moment is to be evaluated # # returns the value of the moment at this sigma moment <- object[[1]] moment.rep <- object[[2]] coefficients.rep <- object[[3]] # evaluate the moment by adding the value at each representation # this is the product of all sigma[i, j]^l[i, j] # if sigma and l are thought of as square matrices and l is the representation moment.value <- 0 for (irep in 1:(dim(moment.rep)[1])) {moment.value <- moment.value + coefficients.rep[irep] * prod(sigma^moment.rep[irep, ])} return(as.vector(moment.value))} print.moment <- function(x, ...){ # function: method to print a moment of the multivariate normal distribution # x is an object of class moment subblank <- function (inputstring) { # remove blanks from a character string outputstring <- sub(" ", "", inputstring) if (outputstring == inputstring) { return(outputstring) } if (outputstring != inputstring) { outputstring <- subblank(outputstring) } return(outputstring) } moment <- x$moment
coef <- as.numeric(x$coefficient) representation <- x$representation
express <- "E["
for (imom in 1:length(moment))
{term <- subblank(paste("X",imom,"^",moment[imom]))
express <- paste(express,term)}
express <- paste(express,"]:")
cat(express, "\n")
print(cbind(coef,representation))

invisible(x)}

make.all.moments <-
function (moment,verbose=TRUE)
{
# make all moments less than or equal to input moment
# put these into the symmoments namespace

subblank <- function (inputstring)
{
# remove blanks from a character string
outputstring <- sub(" ", "", inputstring)
if (outputstring == inputstring) {
return(outputstring)
}
if (outputstring != inputstring) {
outputstring <- subblank(outputstring)
}
return(outputstring)
}

make.moment.name <- function (moment)
{
# returns character name of moment
# moment is the moment as a vector, eg, c(1,2,3)

allchar <- c(0:9,letters,toupper(letters))
ndim <- length(moment)
moment.name <- subblank(paste("m",allchar[moment[1]+1]))
if (ndim > 1)
{for (idim in 2:ndim)
{Ak <- allchar[moment[idim]+1]
moment.name <- subblank(paste(moment.name,Ak))
}
}

return(moment.name)
}

make.moment.vector <- function (moment)
{
# returns character vector of moment
# moment is the moment as a vector, eg, c(1,2,3)

ndim <- length(moment)
moment.vector <- subblank(paste("c(",moment[1]))
if (ndim > 1)
{
moment.vector <- subblank(paste(moment.vector,","))
for (idim in 2:ndim)
{Ak <- moment[idim]
moment.vector <- subblank(paste(moment.vector,Ak))
if (idim < ndim)
{moment.vector <- subblank(paste(moment.vector,","))}
}
}
moment.vector <- subblank(paste(moment.vector,")"))
return(moment.vector)
}
toSquare <- function(L.ut)
{
n <- (-1 + sqrt(1+8*length(L.ut)))/2
L <- 0*diag(n)
start <- 1
for (irow in 1:n)
{
L[irow,] <- c(rep(0,(irow-1)),L.ut[start:(start+n-irow)])
start <- start+n-irow + 1
}
return(L)}

create.envir <- TRUE
if (!exists('symmoments'))
{symmoments <- NULL}
if (exists('symmoments'))
{
if (class(symmoments) == 'environment')
{create.envir <- FALSE}
}
if (create.envir)
{
print('Environment symmoments must exist to receive the moment objects.')
return('Please create it using   symmoments <- new.env()')
}

product = prod(moment+1)
cumproduct = cumprod(moment+1)
mdim = length(moment)
thisone <- rep(0,mdim)
for (mcount in 0:(product-1))
{  #1
remain <- mcount
if (mdim > 1)
{remain <- mcount
for (mindex in mdim:2)
{
denom <- cumproduct[mindex-1]
thisone[mindex] <- trunc(remain/denom)
remain <- remain - denom*thisone[mindex]
}
thisone[1] <- remain
}
if (mdim == 1)
{thisone[1] <- mcount}

if (sum(thisone)%%2 ==0)
{ #2
notmoment <- 0
Tmoment <- make.moment.name(thisone)
if (exists(eval(parse(text=subblank(paste("'",Tmoment,"'")))),envir=symmoments,inherits=FALSE))
{ #3
if (class(eval(parse(text=subblank(paste('symmoments$',Tmoment))))) == 'moment') { notmoment <- 0 if (verbose) {print(paste(subblank(subblank(paste("symmoments$",Tmoment)))," exists"))}
}  # done with this one
if (class(eval(parse(text=subblank(paste("symmoments$",Tmoment))))) != 'moment') {notmoment <- 1} } #3 if (!exists(eval(parse(text=subblank(paste("'",Tmoment,"'")))),envir=symmoments,inherits=FALSE) | notmoment == 1) { #3 sortmoment <- sort(thisone) if (sum(sortmoment==thisone) == length(thisone)) { #4 thisvec <- make.moment.vector(thisone) if (verbose) {print(paste("Starting ",subblank(paste('symmoments$',Tmoment))))}
eval(parse(text=subblank(paste("symmoments$",Tmoment," <- callmultmoments(",thisvec,")")))) } #4 canonical, so just make it if (sum(sortmoment==thisone) != length(thisone)) # unsorted { #4 Smoment <- make.moment.name(sortmoment) if (exists(eval(parse(text=subblank(paste("'",Smoment,"'")))),envir=symmoments,inherits=FALSE)) { #5 if (class(eval(parse(text=subblank(paste("symmoments$",Smoment))))) == 'moment')  # canonical moment exists
{ #6
thisvec <- make.moment.vector(thisone)
if (verbose)
{print(paste("Starting ",subblank(paste("symmoments$",Tmoment))))} eval(parse(text=subblank(paste("symmoments$",Tmoment," <- tounsorted(",thisvec,",symmoments$",Smoment,")")))) } #6 if (class(eval(parse(text=subblank(paste("symmoments$",Smoment))))) != 'moment')
{notmoment <- 1}
} #5

if (!exists(eval(parse(text=subblank(paste("'",Tmoment,"'")))),envir=symmoments,inherits=FALSE) | notmoment == 1)
{ #5
thisvec <- make.moment.vector(sortmoment)
if (verbose)
{print(paste("Starting ",subblank(paste("symmoments$",Smoment))," to create",subblank(paste("symmoments$",Tmoment))))}
eval(parse(text=subblank(paste("symmoments$",Smoment," <- callmultmoments(",thisvec,")")))) thisvec <- make.moment.vector(thisone) if (verbose) {print(paste("Starting ",subblank(paste('symmoments$',Tmoment))))}
eval(parse(text=subblank(paste("symmoments$",Tmoment," <- tounsorted(",thisvec,",symmoments$",Smoment,")"))))
} #5
} #4
} #3
} #2
} #1
return(NULL)
}

toLatex_noncentral <-
function (moment,envir='symmoments')
{
# Compute the Latex representation of a noncentral moment

exists.envir <- FALSE
if (exists(envir))
{
if (class(eval(parse(text=envir))) == 'environment')
{exists.envir <- TRUE}
}
if (!exists.envir)
{return(paste('There is no environment named ', envir)) }

subblank <- function (inputstring)
{
# remove blanks from a character string
outputstring <- sub(" ", "", inputstring)
if (outputstring == inputstring) {
return(outputstring)
}
if (outputstring != inputstring) {
outputstring <- subblank(outputstring)
}
return(outputstring)
}

make.moment.name <- function (moment)
{
# returns character name of moment
# moment is the moment as a vector, eg, c(1,2,3)

allchar <- c(0:9,letters,toupper(letters))
ndim <- length(moment)
moment.name <- subblank(paste("m",allchar[moment[1]+1]))
if (ndim > 1)
{for (idim in 2:ndim)
{Ak <- allchar[moment[idim]+1]
moment.name <- subblank(paste(moment.name,Ak))
}
}

return(moment.name)
}
muproduct.Latex <- function (exponent)
{
# create Latex expression for a mu-product with this exponent
ndim <- length(exponent)
term <- " "
if (exponent[1] > 0)
{
if (exponent[1] == 1)
{term <- "\\mu_{1}"}
if (exponent[1] > 1)
{term <- paste("\\mu_{1}^{",exponent[1],"}")}
}
if (ndim > 1)
{
for (idim in 2:ndim)
{
if (exponent[idim] == 1)
{term <- paste(term,"\\mu_{",idim,"}")}
if (exponent[idim] > 1)
{term <- paste(term,"\\mu_{",idim,"}^{",exponent[idim],"}")}
}
}
if (sum(exponent) == 0)
{term <- " "}
term <- subblank(term)

return(term)
}

doubquote <- subblank(paste("\\", "\\"))
fullmoment <- "E["
for (imoment in(1:(length(moment))) )
{fullmoment <- paste(fullmoment,
"X_{", imoment, "}^{",  moment[imoment],  "}")  }
fullmoment <- paste(fullmoment, "\\mid \\mu, \\Sigma] =", doubquote)
fullmoment <- subblank(fullmoment)

product = prod(moment+1)
cumproduct = cumprod(moment+1)
mdim = length(moment)
thisone <- rep(0,mdim)

# special case: X**0
if (mdim == 1 & sum(moment) == 0)
fullmoment <- paste(fullmoment," 1 \\\\")

if (mdim > 1 | sum(moment) > 0)
{

for (mcount in 0:(product-1))
{
if (mdim > 1)
{remain <- mcount
for (mindex in mdim:2)
{
denom <- cumproduct[mindex-1]
thisone[mindex] <- trunc(remain/denom)
remain <- remain - denom*thisone[mindex]
}
thisone[1] <- remain
}
if (mdim == 1)
{thisone[1] <- mcount}

if (sum(thisone)%%2 ==0)
{
sortmoment <- sort(thisone)
Tmoment <- make.moment.name(sortmoment)
if (eval(parse(text=subblank(paste("!exists('",Tmoment,"',envir=",envir,",inherits=FALSE)")))))
{return(paste(Tmoment,' does not exist in environment ',envir))}

Tmoment <- subblank(paste(envir,'$',Tmoment)) eval(parse(text=paste("thismoment <- ",Tmoment))) if (sum(thismoment$moment == thisone) < length(thismoment$moment)) {thismoment <- tounsorted(thisone,thismoment) } muproduct <- muproduct.Latex(moment-thismoment$moment)

combinproduct <- prod(choose(moment,thismoment$moment)) if (combinproduct == 1) {Acombinproduct <- " "} if (combinproduct > 1) {Acombinproduct <- subblank(paste(combinproduct))} thisLatex <- toLatex(thismoment)[-1] # get rid of initial expression if (length(thisLatex) > 1) {thisLatex[1] <- paste("(",thisLatex[1]) thisLatex[length(thisLatex)] <- subblank(paste(thisLatex[length(thisLatex)],")")) thisLatex <- c(thisLatex," \\\\") } if (length(thisLatex) == 1) {#print(paste(thisLatex,thisLatex=="1")) if (subblank(thisLatex) != "1") {#print(paste("length",length(thisLatex),thisLatex)) thisLatex[1] <- paste("(",thisLatex[1]) thisLatex[length(thisLatex)] <- subblank(paste(thisLatex[length(thisLatex)],")")) thisLatex <- c(thisLatex," \\\\") } else {thisLatex <- c("\\\\")} } if (mcount == 0) {thisLatex <- c(paste(Acombinproduct,"\\;",muproduct),thisLatex)} if (mcount > 0) {thisLatex <- c(paste("+",Acombinproduct,"\\;",muproduct),thisLatex)} fullmoment <-c(fullmoment,thisLatex) } } } return(fullmoment) } evaluate_noncentral <- function (moment,mu,sigma,envir='symmoments') { # Evaluate noncentral moment with mean mu and covariance matrix sigma # moment and mu are vectors # sigma is a vector of the upper diagonal # envir is a character variable containing the name of the environment containing the required central moments exists.envir <- FALSE if (exists(envir)) { if (class(eval(parse(text=envir))) == 'environment') {exists.envir <- TRUE} } if (!exists.envir) {return(paste('There is no environment named ', envir)) } subblank <- function (inputstring) { # remove blanks from a character string outputstring <- sub(" ", "", inputstring) if (outputstring == inputstring) { return(outputstring) } if (outputstring != inputstring) { outputstring <- subblank(outputstring) } return(outputstring) } make.moment.name <- function (moment) { # returns character name of moment # moment is the moment as a vector, eg, c(1,2,3) allchar <- c(0:9,letters,toupper(letters)) ndim <- length(moment) moment.name <- subblank(paste("m",allchar[moment[1]+1])) if (ndim > 1) {for (idim in 2:ndim) {Ak <- allchar[moment[idim]+1] moment.name <- subblank(paste(moment.name,Ak)) } } return(moment.name) } product = prod(moment+1) cumproduct = cumprod(moment+1) mdim = length(moment) thisone <- rep(0,mdim) value <- 0 for (mcount in 0:(product-1)) { if (mdim > 1) {remain <- mcount for (mindex in mdim:2) { denom <- cumproduct[mindex-1] thisone[mindex] <- trunc(remain/denom) remain <- remain - denom*thisone[mindex] } thisone[1] <- remain } if (mdim == 1) {thisone[1] <- mcount} if (sum(thisone)%%2 ==0) { sortmoment <- sort(thisone) Tmoment <- make.moment.name(sortmoment) if (eval(parse(text=subblank(paste("!exists('",Tmoment,"',envir=",envir,",inherits=FALSE)"))))) {return(paste(Tmoment,' does not exist in environment ',envir))} Tmoment <- subblank(paste(envir,'$',Tmoment))
eval(parse(text=paste("thismoment <- ",Tmoment)))
if (mdim > 1)
{
if (sum(thismoment$moment == thisone) < length(thismoment$moment))
{thismoment <- tounsorted(thisone,thismoment)}
}
muproduct <- prod(mu^(moment-thismoment$moment)) combinproduct <- prod(choose(moment,thismoment$moment))
thisvalue <- combinproduct*evaluate.moment(thismoment,sigma)*muproduct
value <- value + thisvalue
}
}
return(value)
}
convert.mpoly <-
function (poly)
{
# convert between a mpoly object and a list giving the corresponding moments and coefficients

if (class(poly)=="mpoly")
{
mpoly.list <- unclass(poly)
matrix.size <- length(mpoly.list)

poly.size <- prod(matrix.size)
coeff <- rep(0,poly.size)
variables <- "coef"
vcount <- 0
for (mcount in 1:poly.size)
{
coeff[mcount] <- unlist(mpoly.list[[mcount]])['coef']
variables <- c(variables,names(unlist(mpoly.list[mcount])))
}
variables <- unique(variables)
if (length(variables) == 1 & variables[1] == 'coef')
{
powers <-  matrix(0,ncol=1,nrow=1)
return(list(coeff=coeff,powers=powers))
}

if (length(variables) > 1)
{
variables <- variables[variables !='coef']
ndim <- length(variables)
powers <- matrix(rep(0,ndim*poly.size),ncol=ndim)

for (mcount in 1:poly.size)
{
thisone <- rep(0,ndim)
thisexp <- unlist(unlist(mpoly.list[mcount]))
thesevars <- c(names(unlist(mpoly.list[mcount])))
nvars <- length(thesevars) - 1
thesevars <- thesevars[thesevars != 'coef']
nvars <- length(thesevars)
if (nvars > 0)
{for (idim in 1:ndim)
{
for (ivar in 1:nvars)
{
if (thesevars[ivar] == variables[idim])
{
thisone[idim] <- thisexp[ivar]
}
}
}
}
powers[mcount,] <- thisone
}
return(list(coeff=coeff,powers=powers))
}
}

if (class(poly)!="mpoly")  # assume conversion from matrix (ie, list) to mpoly
{
n.powers <- dim(poly$powers)[1] if (!is.null(n.powers) | 0==0) { mpoly.list <- vector("list",length=n.powers) moment.size <- dim(poly$powers)[2]

variables <- gsub(" ","",paste("X",1:moment.size))
for (ipower in 1:n.powers)
{
toeval <- " "
thisterm <- poly$powers[ipower,] if (sum(thisterm) == 0) {mpoly.list[1] <- list(c(coef=poly$coeff[ipower]))}

if (sum(thisterm) > 0)
{
thesevars <- variables[thisterm > 0]
thesevalues <- thisterm[thisterm > 0]
toeval <- paste("c(",thesevars[1],"=",thesevalues[1])

if (length(thesevars) > 1)
{
for (ivar in 2:length(thesevars))
{toeval <- paste(toeval,",",thesevars[ivar],"=",thesevalues[ivar])}
}
toeval    <- paste(toeval,", coef = ",poly$coeff[ipower],")") eval(parse(text=paste("mpoly.list[[ipower]] <- ",toeval))) } } output <- mpoly.list return(output) } } } toSquare <- function (L.ut) { n <- (-1 + sqrt(1 + 8 * length(L.ut)))/2 L <- 0 * diag(n) start <- 1 for (irow in 1:n) { L[irow, ] <- c(rep(0, (irow - 1)), L.ut[start:(start + n - irow)]) start <- start + n - irow + 1 } return(L) } convert.multipol <- function (poly) { # convert between a multipol object and a list giving the corresponding moments and coefficients subblank <- function (inputstring) { # remove blanks from a character string outputstring <- sub(" ", "", inputstring) if (outputstring == inputstring) { return(outputstring) } if (outputstring != inputstring) { outputstring <- subblank(outputstring) } return(outputstring) } if (class(poly)=="multipol") { multipol.array <- as.array(poly) array.size <- dim(multipol.array) poly.size <- prod(array.size) moment.size <- length(dim(multipol.array)[dim(multipol.array)>1]) # problem of multipol array component x**0 if (is.null(moment.size) | moment.size==0) {moment.size <- 1} coeff <- rep(0,poly.size) powers <- matrix(nrow=poly.size,ncol=moment.size) cumproduct = cumprod(dim(multipol.array)) thisone <- rep(1,moment.size) nonzero <- 0 for (mcount in 1:poly.size) { thisone <- rep(1,moment.size) remain <- mcount if (moment.size == 1) {thisone <- c(remain)} if (moment.size > 1) { for (mindex in moment.size:2) { thisone[mindex] <- trunc((remain-1)/cumproduct[mindex-1]) remain <- remain - cumproduct[mindex-1]*thisone[mindex] thisone[mindex] <- thisone[mindex] + 1 } } thisone[1] <- remain powers[mcount,] <- thisone - rep(1,moment.size) ### account for indexing in multipol thistemp = paste(thisone[1]) if (moment.size > 1) { for (mindex in 2:moment.size) {thistemp = paste(thistemp,",",thisone[mindex])} } thistemp = subblank(thistemp) coefftemp <- NULL eval(parse(text=paste("coefftemp <- multipol.array[",thistemp,"]"))) if (coefftemp != 0) { nonzero <- nonzero + 1 eval(parse(text=paste("coeff[nonzero] <- multipol.array[",thistemp,"]"))) powers[nonzero,] <- thisone - rep(1,moment.size) ### account for indexing in multipol } } coeff <- coeff[1:nonzero] powers <- as.matrix(powers[1:nonzero,],ncol=dim(powers)[2]) output <- list(coeff=coeff,powers=powers) return(output) } if (class(poly)!="multipol") # assume conversion from matrix (ie, list) to multipol { n.powers <- dim(poly$powers)[1]
moment.size <- dim(poly$powers)[2] if (dim(poly$powers)[1] == 1)
{dimmult <- c(1,poly$powers + rep(1,moment.size))} if (dim(poly$powers)[1]> 1)
{dimmult <- apply(poly$powers, 2, max) + rep(1,moment.size)} n.elements <- prod(dimmult) multipol.array <- array(rep(0,n.elements),dim=dimmult) output <- as.multipol(multipol.array) for (ipower in 1:n.powers) { thistemp = paste(1+poly$powers[ipower,1])
if (moment.size > 1)
{
for (mindex in 2:moment.size)
{thistemp = paste(thistemp,",",(1+poly$powers[ipower,mindex]))} } thistemp = subblank(thistemp) eval(parse(text=paste("multipol.array[",thistemp,"] <- poly$coeff[",ipower,"]")))
}
output <- as.multipol(multipol.array)

return(output)
}

}

evaluate_expected.polynomial <-
function (poly,mu,sigma,envir='symmoments')
{
# compute expected value of a multidimensional polynomial
# poly is either a multipol objects or
# a list with components "powers", which lists the moments involved,
# and "coeff, their coefficients
# mu is the mean as a vector
# sigma is the variance-covariance matrix as the vector of the upper diagonal part, including diagonal
# envir is a character variable containing the name of the environment containing the required central moments

exists.envir <- FALSE
if (exists(envir))
{
if (class(eval(parse(text=envir))) == 'environment')
{exists.envir <- TRUE}
}
if (!exists.envir)
{return(paste('There is no environment named ', envir)) }

temp.poly <- poly
if (class(poly) == "multipol")
{temp.poly <- convert.multipol(poly)}
if (class(poly) == "mpoly")
{temp.poly <- convert.mpoly(poly)}

ndim <- dim(temp.poly$powers)[2] npowers <- dim(temp.poly$powers)[1]
powers <- temp.poly$powers coeff <- temp.poly$coeff

value <- 0
for (imom in 1:npowers)
{
if (coeff[imom] != 0)
{
thisvalue <- coeff[imom]*evaluate_noncentral(powers[imom,],mu,sigma,envir)
value <- value + thisvalue
}
}

return(value)
}

toMoment <-
function (inputobject, tip.label = NULL)
{
if (class(inputobject) != "matching" & !is.matrix(inputobject)) {
if (class(inputobject) != "L-Newick") {
temp <- inputobject
Newick.out <- inputobject
}
if (class(inputobject) == "L-Newick") {
temp <- inputobject$Newick Newick.out <- inputobject$Newick
}
temp <- gsub(";", " ", temp, fixed = TRUE)
temp <- gsub("(", " ", temp, fixed = TRUE)
temp <- gsub(")", " ", temp, fixed = TRUE)
temp <- gsub(",", " ", temp, fixed = TRUE)
temp <- gsub("  ", " ", temp, fixed = TRUE)
tips <- strsplit(temp, " +")[[1]]
tips <- tips[tips != ""]
n <- length(tips)
neworder <- 1:n
if (!is.null(tip.label)) {
for (itips in 1:n) {
neworder[itips] <- grep(tip.label[itips], tips)
}
}
tips <- tips[neworder]
np <- n
temp.tips <- rep(" ", n)
temp.Newick <- gsub(";", " ", inputobject, fixed = TRUE)
for (itips in 1:n) {
temp.tips[itips] <- gsub(" ", "", paste(";", as.character(itips),
";"), fixed = TRUE)
temp.Newick <- gsub(tips[itips], temp.tips[itips],
temp.Newick, fixed = TRUE)
}
n.temp.tips <- n
temp.tip.names <- temp.tips
for (itips in 1:n) {
temp.tips[itips] <- gsub(" ", "", temp.tips[itips],
fixed = TRUE)
}
couples.n <- 1
while (couples.n > 0) {
couples.n <- 0
for (cspec in 1:np) {
for (rspec in 1:np) {
couple <- gsub(" ", "", paste("(", temp.tips[cspec],
",", temp.tips[rspec], ")"), fixed = TRUE)
if (length(grep(couple, temp.Newick)) > 0) {
couples.n <- couples.n + 1
n.temp.tips <- n.temp.tips + 1
temp.tip.names <- c(temp.tip.names, couple)
temp.tips <- c(temp.tips, gsub(" ", "", paste(";",
as.character(n.temp.tips), ";"), fixed = TRUE))
temp.Newick <- gsub(couple, gsub(" ", "",
paste(";", as.character(n.temp.tips), ";"),
fixed = TRUE), temp.Newick, fixed = TRUE)
}
}
}
np <- np + couples.n
}
temp.tip.names.n <- length(temp.tip.names)
L <- 0 * diag(2 * (n - 1))
for (icouple in (n + 1):temp.tip.names.n) {
temp <- gsub("(", " ", temp.tip.names[icouple], fixed = TRUE)
temp <- gsub(")", " ", temp, fixed = TRUE)
temp <- gsub(",", " ", temp, fixed = TRUE)
temp <- gsub(";", " ", temp, fixed = TRUE)
temp <- gsub("  ", " ", temp, fixed = TRUE)
ind <- strsplit(temp, " +")[[1]]
ind <- ind[ind != ""]
indn <- as.integer(ind)
indn <- ind[order(indn)]
irow <- as.integer(indn[1])
icol <- as.integer(indn[2])
L[irow, icol] <- 1
}
L.ut <- NULL
for (irow in 1:(2 * (n - 1))) {
L.ut <- c(L.ut, L[irow, irow:(2 * (n - 1))])
}
SMNM <- list(L = L, L.ut = L.ut, Newick = Newick.out,
tip.label = tips, tip.label.n = n)
class(SMNM) <- "L-matrix"
return(SMNM)
}
if (class(inputobject) == "matching" | is.matrix(inputobject)) {
temp.tip.label <- NULL
if (class(inputobject) == "matching") {
temp.matching <- inputobject$matching[, c(1, 2)] if (!is.null(inputobject$tip.label)) {
temp.tip.label <- inputobject$tip.label } } if (is.matrix(inputobject)) { temp.matching <- inputobject[, c(1, 2)] } nL <- max(temp.matching) n <- 1 + nL/2 L <- 0 * diag(nL) if (is.matrix(temp.matching)) { L[temp.matching] <- 1 } if (!is.matrix(temp.matching)) { L[matrix(temp.matching, nrow = 1)] <- 1 } L.ut <- NULL for (irow in 1:(2 * (n - 1))) { L.ut <- c(L.ut, L[irow, irow:(2 * (n - 1))]) } if (!is.null(temp.tip.label)) { Newick <- toNewick(L, type = "square", tip.label = temp.tip.label) } if (is.null(temp.tip.label)) { Newick <- toNewick(L, type = "square") } if (is.null(tip.label)) { alphab <- gsub(" ", "", paste(rep(letters, times = rep(676, 26)), paste(rep(letters, times = rep(26, 26)), letters))) if (n <= 26) { tips <- letters[1:n] } if (n > 26) { tips <- alphab[1:n] } } if (!is.null(tip.label)) { tips <- tip.label } if (is.null(tip.label) & class(inputobject) == "matching") { if (!is.null(inputobject$tip.label)) {
tips <- inputobject$tip.label } } SMNM <- list(L = L, L.ut = L.ut, Newick = Newick$Newick,
tip.label = tips, tip.label.n = n)
class(SMNM) <- "L-matrix"
return(SMNM)
}
}
toNewick <-
function (L, type = NULL, tip.label = NULL)
{
toSquare <- function (L.ut)
{
n <- (-1 + sqrt(1 + 8 * length(L.ut)))/2
L <- 0 * diag(n)
start <- 1
for (irow in 1:n) {
L[irow, ] <- c(rep(0, (irow - 1)), L.ut[start:(start +
n - irow)])
start <- start + n - irow + 1
}
return(L)
}

if (class(L) == "L-matrix") {
L.matrix <- L$L } if (!is.null(type)) { if (type == "square") { L.matrix <- L } if (type == "ut") { L.matrix <- toSquare(L) } } nL <- dim(L.matrix)[1] n <- 1 + nL/2 alphab <- gsub(" ", "", paste(rep(letters, times = rep(676, 26)), paste(rep(letters, times = rep(26, 26)), letters))) if (!is.null(tip.label)) { tips <- tip.label } if (is.null(tip.label) & class(L) == "L-matrix") { tips <- L$tip.label
}
if (is.null(tip.label) & class(L) != "L-matrix") {
if (n <= 26) {
tips <- letters[1:n]
}
if (n > 26) {
tips <- alphab[1:n]
}
}
nodes <- tips
for (icol in 2:nL) {
irow <- (1:icol)[L.matrix[1:icol, icol] == 1]
if (length(irow) == 1) {
ccol <- nodes[icol]
crow <- nodes[irow]
if (icol <= n) {
couple <- gsub(" ", "", paste("(", crow, ",",
ccol, ")"), fixed = TRUE)
}
if (icol > n) {
couple <- gsub(" ", "", paste("(", ccol, ",",
crow, ")"), fixed = TRUE)
}
nodes <- c(nodes, couple)
}
}
nodes.n <- length(nodes)
temp.Newick <- gsub(" ", "", paste(nodes[nodes.n], ";"))
nodes <- gsub("!", "", nodes, fixed = TRUE)
Newick <- list(Newick = temp.Newick, tip.label = tips, tip.label.n = n,
L = L.matrix, node.label = nodes, nodes.label.n = nodes.n)
class(Newick) <- "L-Newick"
return(Newick)
}
toMatching <-
function (L, type = NULL, tip.label = NULL)
{
toSquare <- function (L.ut)
{
n <- (-1 + sqrt(1 + 8 * length(L.ut)))/2
L <- 0 * diag(n)
start <- 1
for (irow in 1:n) {
L[irow, ] <- c(rep(0, (irow - 1)), L.ut[start:(start +
n - irow)])
start <- start + n - irow + 1
}
return(L)
}

if (class(L) == "L-matrix") {
L.matrix <- L$L if (!is.null(L$tip.label) & is.null(tip.label)) {
tip.label <- L$tip.label } } if (!is.null(type)) { if (type == "square") { L.matrix <- L } if (type == "ut") { L.matrix <- toSquare(L) } } nL <- dim(L.matrix)[1] n <- 1 + nL/2 temp.matching <- matrix(rep(0, 3 * (n - 1)), ncol = 3) if (!is.null(tip.label)) { temp.tip.label <- tip.label } if (is.null(tip.label)) { alphab <- gsub(" ", "", paste(rep(letters, times = rep(676, 26)), paste(rep(letters, times = rep(26, 26)), letters))) if (!is.null(tip.label)) { temp.tip.label <- tip.label } if (is.null(tip.label)) { if (n <= 26) { temp.tip.label <- letters[1:n] } if (n > 26) { temp.tip.label <- alphab[1:n] } } } temp.node.label <- temp.tip.label nrows <- 0 for (icol in 2:nL) { irow <- (1:icol)[L.matrix[1:icol, icol] == 1] if (length(irow) == 1) { ccol <- temp.node.label[icol] crow <- temp.node.label[irow] if (icol <= n) { couple <- gsub(" ", "", paste("(", crow, ",", ccol, ")"), fixed = TRUE) } if (icol > n) { couple <- gsub(" ", "", paste("(", ccol, ",", crow, ")"), fixed = TRUE) } nrows <- nrows + 1 temp.node.label <- c(temp.node.label, couple) temp.matching[nrows, ] <- c(irow, icol, length(temp.node.label)) } } temp.node.label <- gsub("!", "", temp.node.label, fixed = TRUE) match.obj <- list(matching = temp.matching, tip.label = temp.tip.label, node.label = temp.node.label) class(match.obj) <- "matching" return(match.obj) } tounsorted <- function (moment,sorted.moment) { # converts a sorted moment to a specified unsorted moment # eg, m(2,3,5) -> m(5,2,3) # each row is sorted separately # the rows may be ordered differently # moment: noncanonical moment to obtain # moment is in vector form, eg, c(3,1,2) # sorted.moment: canonical moment of class "moment" # this moment is monotone in its powers; toSquare <- function (L.ut) { n <- (-1 + sqrt(1 + 8 * length(L.ut)))/2 L <- 0 * diag(n) start <- 1 for (irow in 1:n) { L[irow, ] <- c(rep(0, (irow - 1)), L.ut[start:(start + n - irow)]) start <- start + n - irow + 1 } return(L) } nestedreps <- function (input.vector, inner.rep, outer.rep) { # replicates input.vector, first by inner.rep, then by outer.rep # input.vector: vector to replicate # inner.rep: count of replicates of elements # outer.rep: count of replicates of resulting vector temp.vector <- NULL for (ivec in 1:length(input.vector)) {temp.vector <- c(temp.vector, rep(input.vector[ivec], inner.rep))} total.vector <- rep(temp.vector, outer.rep) return(total.vector)} output.moment <- sorted.moment output.moment$moment <- moment
lmom <- length(moment)
r <- order(moment)
lrep <- lmom*(lmom+1)/2
representation <- sorted.moment$representation nrep <- dim(representation)[1] overallcoeff <- ((1/2)^(sum(moment)/2))*prod(fact(moment))/fact(sum(moment)/2) m <- matrix(rep(1, lmom^2), nrow=lmom) limits <- c((lmom:1)%*%(m * !(lower.tri(m)))) # sum of row lengths: lmom, lmom + lmom - 1, ... , for use in row_col row_col <- matrix(rep(0, (2 * lrep)), nrow=2) for (icell in (1:lrep)) { row_col[1, icell] <- min((1:lmom)[icell<=limits]) if (row_col[1, icell] == 1){row_col[2, icell] <- icell} if (row_col[1, icell]>1){row_col[2, icell] <- icell - limits[row_col[1, icell] - 1] + row_col[1, icell] - 1 } } # 2x(nm * (nm + 1) / 2 matrix giving rows and columns for each cell # compute here so that they don't have to be calculated each time for (irep in 1:nrep) { utri <- toSquare(representation[irep,]) noncanonical.matrix <- 0*utri for (irow in 1:lmom) { jrow <- r[irow] for (icol in irow:lmom) { jcol <- r[icol] if (jrow <= jcol) {noncanonical.matrix[jrow,jcol] <- utri[irow,icol]} if (jrow > jcol) {noncanonical.matrix[jcol,jrow] <- utri[irow,icol]} } } thisrep <- t(noncanonical.matrix)[t(!lower.tri(noncanonical.matrix))] output.moment$representation[irep,] <- thisrep

#  determine the coefficient for each term based on switching equivalent terms
#          this is taken from callmultmoments

#  "base" gives the number of switches that can be made to each element of the l-matrix
#  diagonal elements are not switchable, but are included to allow subtraction below

base <- c(rep(1, lmom * (lmom + 1) / 2))
base[1] <- 1   # first diagonal element  -  not switchable

totreps <- 1   # total number of transpostions
# if there is only one element, it must be the diagonal, so is not switchable  -  skip
if (lmom > 1){
base[1] <- 1   # first diagonal element
for (cell in 2:length(base))
{
icol = row_col[1, cell]  #  determine if diagonal element
irow = row_col[2, cell]
if (irow == icol){base[cell] <- 1}  # diagonal  -  not switchable
if (icol != irow)
{totreps <- totreps * (1 + thisrep[cell])
base[cell] <- 1 + thisrep[cell]}
}  #  done with computing base and total transpositions (totreps)
}

mcoeff <- 1  #  sum of multinomial coefficients
if (totreps > 1){

#  baserep will represent the lower diagonal (including diagonal)
#  of the augmented matrices
baserep <- matrix(rep(0, totreps * length(base)), nrow=totreps)
basegt1 <- base[base>1]
nbase <- 0
for (ibase in 1:length(base))
{if (base[ibase] > 1)
{nbase = nbase + 1
if (nbase == 1){baserep[, ibase] <- nestedreps(c(0:(basegt1[nbase] - 1)), 1, totreps / prod(basegt1[1:nbase])) }
if (nbase > 1) {baserep[, ibase] <- nestedreps(c(0:(basegt1[nbase] - 1)), prod(basegt1[1:(nbase - 1)]), totreps / prod(basegt1[1:nbase])) }
}
}

#  now go through each transposition
if ( !is.na(totreps) & totreps != 1)
{
mcoeff <- 0
for (jrep in (1:totreps)) # check each transposition
{newrep <- baserep[jrep, ]       #  added lower diagonal elements
addrep <- sort(newrep, decreasing=TRUE)[1:(lmom * (lmom - 1) / 2)]
fulnrep <- c((thisrep - newrep), addrep)
thiscoeff <- ((length(fulnrep))^sum(fulnrep)) * dmultinom(x=fulnrep, prob=rep(1.0, length(fulnrep)))
mcoeff <- mcoeff + thiscoeff
#  the multinomial coefficient is obtained from the multinomial distribution
#  multiply by an appropriate power to get rid of probability
}
}
}
if (is.na(totreps)){mcoeff <- 1}
if (totreps == 1)
{mcoeff <- (length(thisrep))^sum(thisrep) * dmultinom(x=thisrep, prob=rep(1.0, length(thisrep)))}

#  determine full coefficient  -  round because all coefficients should be integers
#                                (Note - this statement has not been proved)

output.moment$coefficients[irep] <- round(overallcoeff * mcoeff) } # end of representations return(output.moment)} integrate.polynomial <- function (poly,mu,sigma,lower=NULL,upper=NULL) { # integrate polynomial moment against MVN # poly: either a multipol objects or # a multipol defined by a list with moment powers and coefficients # mu: mean of multivariate normal as vector # sigma: variance-covariance matrix of multivariate normal # lower, upper: vectors giving limits of integration # if one is NULL, then make it the mean +/- 6 * SD if (is.null(lower)) {lower <- mu - 6*sqrt(diag(sigma))} if (is.null(upper)) {upper <- mu + 6*sqrt(diag(sigma))} thispoly <- poly if (class(poly) == "multipol") {thispoly <- convert.multipol(poly)} if (class(poly) == "mpoly") {thispoly <- convert.mpoly(poly)} ndim <- dim(thispoly$powers)[2]
npowers <- dim(thispoly$powers)[1] powers <- thispoly$powers
coeff <- thispoly$coeff value <- 0 f <- function(x) { y <- x[1]^powers[imom,1] for (idim in 2:ndim) { y <- y*x[idim]^powers[imom,idim] } y <- y*dmvnorm(x,mean=mu,sigma=sigma, log=FALSE) return(y) } for (imom in 1:npowers) {thisvalue <- adaptIntegrate(f,lower,upper)$integral
value <- value + coeff[imom]*thisvalue
}
return(value)}


## Try the symmoments package in your browser

Any scripts or data that you put into this service are public.

symmoments documentation built on May 29, 2017, 10:29 p.m.