if (is.vector(x))  {    n <- 1; x1 <- x[1]; x2 <- x[2] ;x3 <- x[3]; x4 <- x[4];   }  else  {    n <- nrow(x); x1 <- x[,1]; x2 <- x[,2]; x3 <- x[,3]; x4 <- x[,4];  }  d <- 4  sumval <- 0  if (binned)  {    ## drvkde computes include-diagonals estimate    fhatr <- drvkde(x=bin.par$counts, drv=r, bandwidth=sqrt(diag(Sigma)),                       binned=TRUE, range.x=bin.par$range.x, se=FALSE)$est    sumval <- sum(bin.par$counts * n * fhatr)    if (inc == 0)       sumval <- sumval - n*dmvnorm.deriv.4d(x=rep(0,d), r=r, Sigma=Sigma)  }  else  {      for (j in 1:n)    {        y1 <- x1 - x1[j]      y2 <- x2 - x2[j]      y3 <- x3 - x3[j]      y4 <- x4 - x4[j]      sumval <- sumval + sum(dmvnorm.deriv.4d(cbind(y1, y2, y3, y4), Sigma=Sigma, r=r))    }      if (inc==0)      sumval <- sumval - n*dmvnorm.deriv.4d(rep(0,d), Sigma, r)  }    return(sumval)}################################################################################# Double sum  of K(X_i - X_j) used in density derivative estimation - 5-dim### Parameters## x - points to evaluate## Sigma - variance matrix## inc - 0 - exclude diagonals##     - 1 - include diagonals### Returns## Double sum at x###############################################################################dmvnorm.5d.sum <- function(x, Sigma, inc=1){  if (is.vector(x))  {    n <- 1; d <- 4; x1 <- x[1]; x2 <- x[2]; x3 <- x[3]; x4 <- x[4];    x5 <- x[5];   }  else  {    n <- nrow(x); d <- ncol(x); x1 <- x[,1]; x2 <- x[,2]; x3 <- x[,3]; x4 <- x[,4];    x5 <- x[,5];   }    viSigma <- vec(chol2inv(chol(Sigma)))  result <- .C("dmvnorm_5d_sum", as.double(x1), as.double(x2),               as.double(x3), as.double(x4), as.double(x5),               as.double(viSigma), as.double(det(Sigma)), as.integer(n),               as.double(0), PACKAGE="ks")  sumval <- result[[9]]  ## above C function mvnorm_5d_sum only computes the upper triangular half  ## so need to reflect along the diagonal and then subtract appropriate  ## amount to compute whole sum     if (inc == 0)     sumval <- 2*sumval - 2*n*dmvnorm(rep(0,d), rep(0,d), Sigma)  else if (inc == 1)    sumval <- 2*sumval - n*dmvnorm(rep(0,d), rep(0,d), Sigma)     return(sumval)}dmvnorm.deriv.5d.sum <- function(x, Sigma, r, inc=1){  if (is.vector(x))  {    n <- 1; d <- 4; x1 <- x[1]; x2 <- x[2]; x3 <- x[3]; x4 <- x[4];    x5 <- x[5];   }  else  {    n <- nrow(x); d <- ncol(x); x1 <- x[,1]; x2 <- x[,2]; x3 <- x[,3]; x4 <- x[,4];    x5 <- x[,5];   }  d <- 5  sumval <- 0  for (j in 1:n)  {      y1 <- x1 - x1[j]    y2 <- x2 - x2[j]    y3 <- x3 - x3[j]    y4 <- x4 - x4[j]    y5 <- x5 - x5[j]    sumval <- sumval + sum(dmvnorm.deriv.5d(cbind(y1, y2, y3, y4, y5),Sigma,r))  }    if (inc==0)    sumval <- sumval - n*dmvnorm.deriv.5d(rep(0,d), Sigma, r)      return(sumval)}################################################################################# Double sum  of K(X_i - X_j) used in density derivative estimation - 6-dim### Parameters## x - points to evaluate## Sigma - variance matrix## inc - 0 - exclude diagonals##     - 1 - include diagonals### Returns## Double sum at x###############################################################################dmvnorm.6d.sum <- function(x, Sigma, inc=1){  if (is.vector(x))  {    n <- 1; d <- 4; x1 <- x[1]; x2 <- x[2]; x3 <- x[3]; x4 <- x[4];    x5 <- x[5]; x6 <- x[6];   }  else  {    n <- nrow(x); d <- ncol(x); x1 <- x[,1]; x2 <- x[,2]; x3 <- x[,3]; x4 <- x[,4];    x5 <- x[,5]; x6 <- x[,6];    }    viSigma <- vec(chol2inv(chol(Sigma)))  result <- .C("dmvnorm_6d_sum", as.double(x1), as.double(x2),               as.double(x3), as.double(x4), as.double(x5), as.double(x6),               as.double(viSigma), as.double(det(Sigma)), as.integer(n),               as.double(0), PACKAGE="ks")  sumval <- result[[10]]  ## above C function mvnorm_6d_sum only computes the upper triangular half  ## so need to reflect along the diagonal and then subtract appropriate  ## amount to compute whole sum     if (inc == 0)     sumval <- 2*sumval - 2*n*dmvnorm(rep(0,d), rep(0,d), Sigma)  else if (inc == 1)    sumval <- 2*sumval - n*dmvnorm(rep(0,d), rep(0,d), Sigma)     return(sumval)}dmvnorm.deriv.6d.sum <- function(x, Sigma, r, inc=1){  if (is.vector(x))  {    n <- 1; d <- 4; x1 <- x[1]; x2 <- x[2]; x3 <- x[3]; x4 <- x[4];    x5 <- x[5]; x6 <- x[6];   }  else  {    n <- nrow(x); d <- ncol(x); x1 <- x[,1]; x2 <- x[,2]; x3 <- x[,3]; x4 <- x[,4];    x5 <- x[,5]; x6 <- x[,6];  }  d <- 6  sumval <- 0  for (j in 1:n)  {      y1 <- x1 - x1[j]    y2 <- x2 - x2[j]    y3 <- x3 - x3[j]    y4 <- x4 - x4[j]    y5 <- x5 - x5[j]    y6 <- x6 - x6[j]    sumval <- sumval + sum(dmvnorm.deriv.6d(cbind(y1, y2, y3, y4, y5, y6),Sigma,r))  }    if (inc==0)    sumval <- sumval - n*dmvnorm.deriv.6d(rep(0,d), Sigma, r)      return(sumval)}################################################################################# Compute moments of multivariate normal mixture###############################################################################moments.mixt <- function (mus, Sigmas, props){  if (!(identical(all.equal(sum(props), 1), TRUE)))    stop("Proportions don't sum to one\n")  d <- ncol(Sigmas)  k <- length(props)  mn <- rep(0, d)  va <- matrix(0, nr=d, nc=d)  for (i in 1:k)  {    mn <- mn + props[i] * mus[i,]    va <- va + props[i] * (Sigmas[((i-1)*d+1):(i*d),] + mus[i,] %*% t(mus[i,]))  }   va <- va + mn %*% t(mn)  return( list(mean=mn, var=va))}################################################################################# Creates plots of mixture density functions### Parameters## mus - means## Sigmas - variances## props - vector of proportions of each mixture component ## dfs - degrees of freedom## dist - "normal" - normal mixture##      - "t" - t mixture## ...###############################################################################plotmixt <- function(mus, Sigmas, props, dfs, dist="normal", ...){  if (ncol(Sigmas)==2)    plotmixt.2d(mus=mus, Sigmas=Sigmas, props=props, dfs=dfs, dist=dist, ...)  else if (ncol(Sigmas)==3)    plotmixt.3d(mus=mus, Sigmas=Sigmas, props=props, dfs=dfs, dist=dist, ...) }plotmixt.2d <- function(mus, Sigmas, props, dfs, dist="normal",    xlim, ylim, gridsize, display="slice", cont=c(25,50,75), abs.cont,    lty, xlab="x", ylab="y", zlab="Density function",    theta=-30, phi=40, d=4, add=FALSE, drawlabels=TRUE, nrand=1e5, ...){  dist <- tolower(substr(dist,1,1))  maxSigmas <- 4*max(Sigmas)  if (is.vector(mus))    mus <- as.matrix(t(mus))  if (missing(xlim))    xlim <- c(min(mus[,1]) - maxSigmas, max(mus[,1]) + maxSigmas)  if (missing(ylim))    ylim <- c(min(mus[,2]) - maxSigmas, max(mus[,2]) + maxSigmas)  if (missing(gridsize))    gridsize <- rep(51,2)                x <- seq(xlim[1], xlim[2], length=gridsize[1])  y <- seq(ylim[1], ylim[2], length=gridsize[2])  xy <- permute(list(x, y))  d <- ncol(Sigmas)    if (dist=="n")    dens <- dmvnorm.mixt(xy, mu=mus, Sigma=Sigmas, props=props)    else if (dist=="t")    dens <- dmvt.mixt(xy, mu=mus, Sigma=Sigmas, props=props, dfs=dfs)  dens.mat <- matrix(dens, nc=length(x), byrow=FALSE)     disp <- substr(display,1,1)  if (disp=="p")    persp(x, y, dens.mat, theta=theta, phi=phi, d=d, xlab=xlab, ylab=ylab,          zlab=zlab, ...)  else if (disp=="s")  {    if (dist=="n")    {      x.rand <- rmvnorm.mixt(n=nrand, mus=mus, Sigmas=Sigmas, props=props)      dens.rand <- dmvnorm.mixt(x.rand, mus=mus, Sigmas=Sigmas, props=props)    }    else if (dist=="t")    {      x.rand <- rmvt.mixt(n=nrand, mus=mus, Sigmas=Sigmas, props=props, dfs=dfs)      dens.rand <- dmvt.mixt(x.rand, mus=mus, Sigmas=Sigmas, props=props, dfs=dfs)    }        if (missing(lty))      lty <- 1    if (missing(abs.cont))      hts <- quantile(dens.rand, prob=(100 - cont)/100)    else      hts <- abs.cont        if (!add)      plot(x, y, type="n", xlab=xlab, ylab=ylab, xlim=xlim, ylim=ylim, ...)            for (i in 1:length(hts))     {      scale <- cont[i]/hts[i]      if (missing(abs.cont))        contour(x, y, dens.mat*scale, level=hts[i]*scale, add=TRUE,                drawlabels=drawlabels, lty=lty, ...)      else        contour(x, y, dens.mat, level=hts[i], add=TRUE,                drawlabels=drawlabels, lty=lty, ...)    }  }    else if (disp=="i")    image(x, y, dens.mat, xlab=xlab, ylab=ylab, ...)  else if (disp=="f")    filled.contour(x, y, dens.mat, xlab=xlab, ylab=ylab, ...)    }plotmixt.3d <- function(mus, Sigmas, props, dfs, cont=c(25,50,75), abs.cont,    dist="normal", xlim, ylim, zlim, gridsize, alphavec, colors, add=FALSE, nrand=1e5, ...){  require(rgl)  require(misc3d)  d <- 3  dist <- tolower(substr(dist,1,1))  maxSigmas <- 3.7*max(Sigmas)  if (is.vector(mus))    mus <- as.matrix(t(mus))    if (missing(xlim))    xlim <- c(min(mus[,1]) - maxSigmas, max(mus[,1]) + maxSigmas)  if (missing(ylim))    ylim <- c(min(mus[,2]) - maxSigmas, max(mus[,2]) + maxSigmas)  if (missing(zlim))    zlim <- c(min(mus[,3]) - maxSigmas, max(mus[,3]) + maxSigmas)    if (missing(gridsize))    gridsize <- rep(51,d)    x <- seq(xlim[1], xlim[2], length=gridsize[1])  y <- seq(ylim[1], ylim[2], length=gridsize[2])  z <- seq(zlim[1], zlim[2], length=gridsize[3])  xy <- permute(list(x,y))  dens.array <- array(0, dim=gridsize)    for (i in 1:length(z))  {    if (dist=="n")      dens <- dmvnorm.mixt(cbind(xy, z[i]), mu=mus, Sigma=Sigmas, props=props)    else if (dist=="t")      dens <- dmvt.mixt(cbind(xy, z[i]), mu=mus, Sigma=Sigmas, dfs=dfs, props=props)        dens.mat <- matrix(dens, nc=length(x), byrow=FALSE)    dens.array[,,i] <- dens.mat  }    if (dist=="n")  {      x.rand <- rmvnorm.mixt(n=nrand, mus=mus, Sigmas=Sigmas, props=props)    dens.rand <- dmvnorm.mixt(x.rand, mus=mus, Sigmas=Sigmas, props=props)  }  else if (dist=="t")  {    x.rand <- rmvt.mixt(n=nrand, mus=mus, Sigmas=Sigmas, props=props, dfs=dfs)    dens.rand <- dmvt.mixt(x.rand, mus=mus, Sigmas=Sigmas, props=props, dfs=dfs)  }  if (missing(abs.cont))    hts <- quantile(dens.rand, prob = (100 - cont)/100)  else    hts <- abs.cont  nc <- length(hts)  if (missing(colors))    colors <- rev(heat.colors(nc))    if (missing(alphavec))    alphavec <- seq(0.1,0.5,length=nc)  plot3d(x, y, z, type="n", add=add, ...)  for (i in 1:nc)   {    ##scale <- cont[i]/hts[i]    contour3d(dens.array, level=hts[nc-i+1],x, y, z, add=TRUE, color=colors[i],             alpha=alphavec[i], ...)  }}################################################################################# Multivariate t - density values### Parameters## x - points to compute density     ## mu - vector of means ## Sigma - dispersion matrix## df - degrees of freedom### Returns## Value of multivariate t density at x###############################################################################dmvt <- function(x, mu, Sigma, df){     if(is.vector(x))    x <- matrix(x, ncol=length(x))  d <- ncol(Sigma)  detSigma <- det(Sigma)   dens <- (1+ mahalanobis(x, center=mu, cov=Sigma)/df)^(-(d+df)/2)  dens <- dens * gamma((df+d)/2) / ((df*pi)^(d/2)*gamma(df/2)*detSigma^(1/2))    return(dens)}################################################################################# Multivariate t mixture - density values### Parameters## x - points to compute density at    ## mus - vector of means ## Sigmas - dispersion matrices## dfs - degrees of freedom## props - vector of mixing proportions### Returns## Value of multivariate t mixture density at x###############################################################################dmvt.mixt <- function(x, mus, Sigmas, dfs, props){  if (!(identical(all.equal(sum(props), 1), TRUE)))    stop("Proportions don't sum to one\n")  else if (length(dfs) != length(props))    stop("Length of df and mixing proportions vectors not equal")    ## single component mixture  if (identical(all.equal(props[1], 1), TRUE))    dens <- dmvt(x, mu=mus, Sigma=Sigmas, df=dfs)    ## multiple component mixture  else     {       if (is.vector(mus)) d <- length(mus)    else d <- ncol(mus)    k <- length(props)    dens <- 0          for (i in 1:k)      dens <- dens+props[i]*dmvt(x,mu=mus[i,],Sigma=Sigmas[((i-1)*d+1):(i*d),],                                 df=dfs[i])  }    return(dens)}################################################################################# Multivariate t mixture - random sample## ## Parameters## n - number of samples## mus - means ## Sigmas - matrix of dispersion matrices## dfs - vector of degrees of freedom## props - vector of mixing proportions ## ## Returns## Vector of n observations from the t mixture###############################################################################rmvt.mixt <- function(n=100, mus=c(0,0), Sigmas=diag(2), dfs=7, props=1){  if (!(identical(all.equal(sum(props), 1), TRUE)))      stop("Proportions don't sum to one\n")  else if (length(dfs) != length(props))    stop("Length of df and mixing proportions vectors not equal")    ## single component mixture  if (identical(all.equal(props[1], 1), TRUE))  {    rand <- rmvt(n=n, sigma=Sigmas, df=dfs)    for (i in 1:length(mus))      rand[,i] <- rand[,i] + mus[i]  }    ## multiple component mixture  else  {    k <- length(props)    d <- ncol(Sigmas)    n.samp <- sample(1:k, n, replace=TRUE, prob=props)     n.prop <- numeric(0)    ## compute number to be drawn from each component     for (i in 1:k)      n.prop <- c(n.prop, sum(n.samp == i))    ## generate random samples from each component    rand <- numeric(0)      for (i in 1:k)    {      if (n.prop[i] > 0)      {          rand.temp<-rmvt(n=n.prop[i],sigma=Sigmas[((i-1)*d+1):(i*d),],df=dfs[k])        for (j in 1:length(mus[k,]))          rand.temp[,j] <- rand.temp[,j] + mus[i,j]               rand <- rbind(rand, rand.temp)      }    }  }    return(rand[sample(n),])}