################################################################################# Univariate mixture normal densities###############################################################################rnorm.mixt <- function(n=100, mus=0, sigmas=1, props=1, mixt.label=FALSE){  if (!(identical(all.equal(sum(props), 1), TRUE)))    stop("Proportions don't sum to one\n")  ### single component mixture  if (identical(all.equal(props[1], 1), TRUE))  {    if (mixt.label)      rand <- cbind(rnorm(n=n, mean=mus, sd=sigmas), rep(1, n))    else      rand <- rnorm(n=n, mean=mus, sd=sigmas)  }  ### multiple component mixture  else  {    k <- length(props)    n.samp <- sample(1:k, n, replace=TRUE, prob=props)     n.prop <- numeric(0)    ## alternative method for component membership    ##runif.memb <- runif(n=n)    ##memb <- findInterval(runif.memb, c(0,cumsum(props)), rightmost.closed=TRUE)    ##n.prop <- table(memb)        ## compute number taken from each mixture    for (i in 1:k)      n.prop <- c(n.prop, sum(n.samp == i))        rand <- numeric(0)        for (i in 1:k) ##for (i in as.numeric(rownames(n.prop)))    {      ## compute random sample from normal mixture component      if (n.prop[i] > 0)        if (mixt.label)          rand <- rbind(rand, cbind(rnorm(n=n.prop[i], mean=mus[i], sd=sigmas[i]), rep(i, n.prop[i])))        else          rand <- c(rand, rnorm(n=n.prop[i], mean=mus[i], sd=sigmas[i]))    }  }  if (mixt.label)    return(rand[sample(n),])  else    return(rand[sample(n)])}dnorm.mixt <- function(x, mus=0, sigmas=1, props=1){  if (!(identical(all.equal(sum(props), 1), TRUE)))    stop("Proportions don't sum to one\n")  ## single component mixture  if (identical(all.equal(props[1], 1), TRUE))    dens <- dnorm(x, mean=mus, sd=sigmas)  ## multiple component mixture  else     {       k <- length(props)    dens <- 0    ## sum of each normal density value from each component at x      for (i in 1:k)      dens <- dens + props[i]*dnorm(x, mean=mus[i], sd=sigmas[i])  }    return(dens)}   ################################################################################# Partial derivatives of the univariate normal (mean 0) ## ## Parameters## x - points to evaluate at## sigma - std deviation## r - derivative index ### Returns## r-th derivative at x###############################################################################dnorm.deriv <- function(x, mu=0, sigma=1, r=0){  phi <- dnorm(x, mean=mu, sd=sigma)   x <- (x - mu)    if (r==0)    return(phi)  else if (r==1)    derivt <- -x/sigma^2*phi  else if (r==2)    derivt <- (x^2-sigma^2)/sigma^4*phi  else if (r==3)    derivt <- -(x^3 - 3*x*sigma^2)/sigma^6*phi  else if (r==4)    derivt <- (x^4 - 6*x^2*sigma^2 + 3*sigma^4)/sigma^8*phi  else if (r==5)    derivt <- -(x^5 - 10*x^3*sigma^2 + 15*x*sigma^4)/sigma^10*phi  else if (r==6)    derivt <- (x^6 - 15*x^4*sigma^2 + 45*x^2*sigma^4 - 15*sigma^6)/sigma^12*phi  else if (r==7)    derivt <- -(x^7 - 21*x^5*sigma^2 + 105*x^3*sigma^4 - 105*x*sigma^6)/sigma^14*phi  else if (r==8)    derivt <- (x^8 - 28*x^6*sigma^2 + 210*x^4*sigma^4 - 420*x^2*sigma^6 + 105*sigma^8)/sigma^16*phi  else if (r==9)    derivt <- -(x^9 - 36*x^7*sigma^2 + 378*x^5*sigma^4 - 1260*x^3*sigma^6 + 945*x*sigma^8)/sigma^18*phi  else if (r==10)    derivt <- (x^10 - 45*x^8*sigma^2 + 630*x^6*sigma^4 - 3150*x^4*sigma^6 + 4725*x^2*sigma^8 - 945*sigma^10)/sigma^20*phi    if (r > 10)    stop ("Up to 10th order derivatives only")      return(derivt)}################################################################################# Double sum  of K(X_i - X_j) used in density derivative estimation### Parameters## x - points to evaluate## Sigma - variance matrix## inc - 0 - exclude diagonals##     - 1 - include diagonals### Returns## Double sum at x###############################################################################dnorm.sum <- function(x, sigma=1, inc=1, binned=FALSE, bin.par){  d <- 1  if (binned)  {    if (missing(bin.par)) bin.par <- binning(x, h=sigma)      n <- sum(bin.par$counts)    fhatr <- kde.binned(bin.par=bin.par, h=sigma)    sumval <- sum(bin.par$counts * n * fhatr$estimate)    ##sumval <- bkfe(x=bin.par$counts, bandwidth=sigma, drv=0, binned=TRUE, range.x=range(bin.par$eval.points))    if (inc == 0)       sumval <- sumval - n*dnorm.deriv(x=0, mu=0, r=0, sigma=sigma)  }  else  {    n <- length(x)    sumval <- 0    for (i in 1:n)      sumval <- sumval + sum(dnorm(x=x[i] - x, mean=0, sd=sigma))        if (inc == 0)       sumval <- sumval - n*dnorm(x=0, mean=0, sd=sigma)   }     return(sumval)}dnorm.deriv.sum <- function(x, sigma, r, inc=1, binned=FALSE, bin.par, kfe=FALSE){  if (binned)  {    if (missing(bin.par)) bin.par <- binning(x, h=sigma, supp=4+r)      fhatr <- kdde.binned(bin.par=bin.par, h=sigma, r=r)    n <- sum(bin.par$counts)    sumval <- sum(bin.par$counts * n * fhatr$estimate)    ##sumval <- n*bkfe(x=bin.par$counts, bandwidth=sigma, drv=r, binned=TRUE, range.x=range(bin.par$eval.points))  }  else  {    n <- length(x)    sumval <- 0    for (i in 1:n)      sumval <- sumval + sum(dnorm.deriv(x=x[i] - x, mu=0, sigma=sigma, r=r))   }  if (inc == 0)     sumval <- sumval - n*dnorm.deriv(x=0, mu=0, sigma=sigma, r=r)  if (kfe)    if (inc==1) sumval <- sumval/n^2    else sumval <- sumval/(n*(n-1))    return(sumval)  }################################################################################ Multivariate normal densities and derivatives################################################################################################################################################################ Multivariate normal mixture - random sample## ## Parameters## n - number of samples## mus - matrix of means (each row is a vector of means from each component##       density)## Sigmas - matrix of covariance matrices (every d rows is a covariance matrix ##          from each component density) ## props - vector of mixing proportions ## ## Returns## Vector of n observations from the normal mixture ###############################################################################rmvnorm.mixt <- function(n=100, mus=c(0,0), Sigmas=diag(2), props=1, mixt.label=FALSE){  if (!(identical(all.equal(sum(props), 1), TRUE)))    stop("Proportions don't sum to one\n")    #if (is.vector(Sigmas))  ##  return(rnorm.mixt(n=n, mus=mus, sigmas=Sigmas, props=props))    ### single component mixture  if (identical(all.equal(props[1], 1), TRUE))   if (mixt.label)     rand <- cbind(rmvnorm(n=n, mean=mus, sigma=Sigmas), rep(1, n))   else     rand <- cbind(rmvnorm(n=n, mean=mus, sigma=Sigmas))      ### 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 taken from each mixture    for (i in 1:k)      n.prop <- c(n.prop, sum(n.samp == i))        rand <- numeric(0)        for (i in 1:k)    {      ## compute random sample from normal mixture component      if (n.prop[i] > 0)      {               if (mixt.label)          rand <- rbind(rand, cbind(rmvnorm(n=n.prop[i], mean=mus[i,], sigma=Sigmas[((i-1)*d+1) : (i*d),]), rep(i, n.prop[i])))        else          rand <- rbind(rand, rmvnorm(n=n.prop[i], mean=mus[i,], sigma=Sigmas[((i-1)*d+1) : (i*d),]))      }        }  }  return(rand[sample(n),])}################################################################################# Multivariate normal mixture - density values## ## Parameters## x - points to compute density at ## mus - matrix of means## Sigmas - matrix of covariance matrices ## props - vector of mixing proportions ## ## Returns## Density values from the normal mixture (at x)###############################################################################dmvnorm.mixt <- function(x, mus, Sigmas, props=1){    if (!(identical(all.equal(sum(props), 1), TRUE)))    stop("Proportions don't sum to one\n")  if (is.vector(x)) d <- length(x)  else d <- ncol(x)    if (missing(mus)) mus <- rep(0,d)  if (missing(Sigmas)) Sigmas <- diag(d)  ##if (missing(deriv)) deriv <- rep(0,d)     ## single component mixture  if (identical(all.equal(props[1], 1), TRUE))    dens <- dmvnorm(x=x, mean=mus, sigma=Sigmas)  ## multiple component mixture  else     {       k <- length(props)    dens <- 0    ## sum of each normal density value from each component at x      for (i in 1:k)      dens <- dens + props[i]*dmvnorm(x, mean=mus[i,], sigma=Sigmas[((i-1)*d+1):(i*d),])  }    return(dens)}   ################################################################################# Partial derivatives of the multivariate normal## ## Parameters## x - points to evaluate at## Sigma - variance## r - derivative index ### Returns## r-th derivative at x################################################################################## for diagonal Sigmadmvnorm.deriv.diag <- function(x, mu, Sigma, r){  if (is.vector(x))    x <- t(as.matrix(x))  if (is.data.frame(x)) x <- as.matrix(x)  d <- ncol(x)  n <- nrow(x)    if (missing(mu)) mu <- rep(0,d)  if (missing(Sigma)) Sigma <- diag(d)  for (i in 1:n)    x[i,] <- x[i,] - mu      if (d==2) return(dmvnorm.deriv.2d(x=x, Sigma=Sigma, r=r))  if (d==3) return(dmvnorm.deriv.3d(x=x, Sigma=Sigma, r=r))  if (d==4) return(dmvnorm.deriv.4d(x=x, Sigma=Sigma, r=r))  if (d==5) return(dmvnorm.deriv.5d(x=x, Sigma=Sigma, r=r))  if (d==6) return(dmvnorm.deriv.6d(x=x, Sigma=Sigma, r=r))}### for general Sigmadmvnorm.deriv <- function(x, mu, Sigma, r, Sdr.mat){     if (is.vector(x))    x <- t(as.matrix(x))  if (is.data.frame(x)) x <- as.matrix(x)  d <- ncol(x)  n <- nrow(x)    sumr <- sum(r)    if (missing(mu)) mu <- rep(0,d)  if (missing(Sigma)) Sigma <- diag(d)  for (i in 1:n)      x[i,] <- x[i,] - mu    if (missing(Sdr.mat) & d >=2)    Sdr.mat <- Sdr(d=d, r=sumr)  dens <- dmvnorm(x=x, mean=mu, sigma=Sigma)   vSigma <- vec(Sigma)  Sigmainv <- chol2inv(chol(Sigma))  mvh <- matrix(0, nrow=n, ncol=d^sumr)   if (sumr==0)    mvh <- dens  if (sumr==1)    mvh <- x*dens  if (sumr==2)  {    Sinv <- (Sigmainv %x% Sigmainv) %*%  Sdr.mat    for (i in 1:n)      mvh[i,]  <-  Sinv %*% ((x[i,] %x% x[i,]) - vSigma) * dens[i]    ind.mat <- K.sum(diag(d), diag(d))  }  if (sumr==3)  {    Sinv <- K.pow(Sigmainv,3) %*% Sdr.mat    for (i in 1:n)      mvh[i,] <- Sinv %*% (K.pow(x[i,], 3) - 3*x[i,] %x% vSigma) * dens[i]    ind.mat <- K.sum(diag(d), K.sum(diag(d), diag(d)))  }  if (sumr==4)  {        Sinv <- K.pow(Sigmainv,4) %*% Sdr.mat    for (i in 1:n)      mvh[i,] <- Sinv %*% (K.pow(x[i,], 4) - 6*K.pow(x[i,],2) %x% vSigma + 3*vSigma %x% vSigma) * dens[i]    ind.mat <- K.sum(diag(d), K.sum(diag(d), K.sum(diag(d), diag(d))))  }  if (sumr==5)  {    Sinv <- K.pow(Sigmainv,5) %*% Sdr.mat    for (i in 1:n)      mvh[i,] <- Sinv %*% (K.pow(x[i,], 5) - 10*K.pow(x[i,],3) %x% vSigma + 15*x[i,] %x% K.pow(vSigma,2)) * dens[i]    ind.mat <- K.sum(diag(d),K.sum(diag(d), K.sum(diag(d), K.sum(diag(d), diag(d)))))  }  if (sumr==6)  {    Sinv <- K.pow(Sigmainv,6) %*% Sdr.mat    for (i in 1:n)      mvh[i,] <- Sinv %*% (K.pow(x[i,],6) - 15*K.pow(x[i,],4) %x% vSigma + 45*K.pow(x[i,],2) %x% K.pow(vSigma,2) - 15*K.pow(vSigma,3)) * dens[i]    ind.mat <- K.sum(diag(d),K.sum(diag(d), K.sum(diag(d), K.sum(diag(d), K.sum(diag(d), diag(d))))))  }  if (sumr==7)  {    Sinv <- K.pow(Sigmainv,7) %*% Sdr.mat    for (i in 1:n)      mvh[i,] <- Sinv %*% (K.pow(x[i,],7) - 21*K.pow(x[i,],5) %x% vSigma + 105*K.pow(x[i,],3) %x% K.pow(vSigma,2)) * dens[i]    ind.mat <- K.sum(diag(d), K.sum(diag(d),K.sum(diag(d), K.sum(diag(d), K.sum(diag(d), K.sum(diag(d), diag(d)))))))  }  if (sumr==8)  {    Sinv <- K.pow(Sigmainv,8) %*% Sdr.mat    for (i in 1:n)      mvh[i,] <- Sinv %*% (K.pow(x[i,],8) - 28*K.pow(x[i,],6) %x% vSigma + 210*K.pow(x[i,],4) %x% K.pow(vSigma,2) - 420*K.pow(x[i,],2) %x% K.pow(vSigma,3) + 105*K.pow(vSigma, 4)) * dens[i]    ind.mat <- K.sum(diag(d), K.sum(diag(d), K.sum(diag(d),K.sum(diag(d), K.sum(diag(d), K.sum(diag(d), K.sum(diag(d), diag(d))))))))  }  if (sumr > 8)    stop ("Up to 8th order derivatives only")  mvh <- (-1)^sumr*mvh[,!duplicated(ind.mat)]  deriv.ind <- unique(ind.mat)  if (length(r)>1)  {    which.deriv <- which.mat(r, deriv.ind)    if (is.vector(mvh)) return(mvh[which.deriv])    else return(mvh[,which.deriv])  }  else    return(list(deriv=mvh, deriv.ind=deriv.ind))}################################################################################# Double sum  of K(X_i - X_j) used in density derivative estimation### Parameters## x - points to evaluate## Sigma - variance matrix## inc - 0 - exclude diagonals##     - 1 - include diagonals### Returns## Double sum at x##############################################################################dmvnorm.sum <- function(x, Sigma, inc=1, binned=FALSE, bin.par, diff=FALSE){  if (binned)  {    if (!is.diagonal(Sigma))      stop("Binned estimation defined for diagonal Sigma only")    if (missing(bin.par)) bin.par <- binning(x, H=Sigma)          fhatr <- kde.binned(bin.par=bin.par, H=Sigma)$estimate     n <- sum(bin.par$counts)    d <- ncol(Sigma)    sumval <- sum(bin.par$counts * n * fhatr)    if (inc == 0)       sumval <- sumval - n*dmvnorm(x=rep(0,d), mean=rep(0,d), sigma=Sigma)  }  else  {    ### Need to rewrite this for d <=6    d <- ncol(Sigma)    if (d==2) sumval <- dmvnorm.2d.sum(x=x, Sigma=Sigma, inc=inc)    if (d==3) sumval <- dmvnorm.3d.sum(x=x, Sigma=Sigma, inc=inc)    if (d==4) sumval <- dmvnorm.4d.sum(x=x, Sigma=Sigma, inc=inc)    if (d==5) sumval <- dmvnorm.5d.sum(x=x, Sigma=Sigma, inc=inc)    if (d==6) sumval <- dmvnorm.6d.sum(x=x, Sigma=Sigma, inc=inc)    if (d>6)    {      if(!diff)      {        n <- nrow(x)        d <- ncol(x)        difs <- differences(x, upper=TRUE)      }      else      {        n <- (-1 + sqrt(1+8*nrow(x)))/2        d <- ncol(x)        difs <- x      }      sumval <- sum(dmvnorm(difs, mean=rep(0,d), sigma=Sigma))            if (inc==0)        sumval <- 2*sumval - 2*n*dmvnorm(rep(0,d), mean=rep(0,d), sigma=Sigma)      if (inc==1)        sumval <- 2*sumval - n*dmvnorm(rep(0,d), mean=rep(0,d), sigma=Sigma)     }    }  return(sumval)}dmvnorm.deriv.sum <- function(x, Sigma, r, inc=1, binned=FALSE, bin.par, diff=FALSE, kfe=FALSE){  if (binned)  {    if (!is.diagonal(Sigma))      stop("Binned estimation defined for diagonal Sigma only")    if (missing(bin.par)) bin.par <- binning(x, H=Sigma)      fhatr <- kdde.binned(bin.par=bin.par, H=Sigma, r=r)$estimate     n <- sum(bin.par$counts)    d <- ncol(Sigma)    sumval <- sum(bin.par$counts * n * fhatr)  }  else  {    if(!diff)    {      n <- nrow(x)