################################################################################ Exact MISE for normal mixtures################################################################################# nu, gamma.r, gamma.r2 written by Jose Chacon 10/2008nu <- function(r,A){ ###Using the recursive formula provided in Kan (2008)  A <- solve(A)  ei <- eigen(A)$values  tr <- numeric(r)  for(p in 1:r)    tr[p] <- sum(ei^p)  nu <- 1  if (r>=1)  {    for(p in 1:r)    {      a<-sum(tr[1:p]*rev(nu))/(2*p)      nu<-c(nu,a)    }  }   return(factorial(r)*2^r*nu[r+1])   }   ## gamma functional for MISE gamma.r <- function(mu, Sigma, d, r, Sd2r){  Sigmainv <- chol2inv(chol(Sigma))  w <- vec(K.pow(Sigmainv %*% Sigmainv, r)) %*% Sd2r  v <- rep(0,length=d^(2*r))  for(j in 0:r)    v <- v + ((-1)^j*OF(2*j)*choose(2*r, 2*j))*(K.pow(mu,2*r-2*j)%x%K.pow(vec(Sigma),j))  gamr <- (-1)^r*dmvnorm(mu,mean=rep(0,d),sigma=Sigma)*sum(w %*% v)    return(gamr)}## gamma functional for AMISE gamma.r2 <- function(mu, Sigma, d, r, Sd2r4, H){    Sigmainv <- chol2inv(chol(Sigma))  w <- vec(K.pow(Sigmainv %*% Sigmainv, r)) %x% vec(K.pow(Sigmainv %*% H %*% Sigmainv, 2)) %*% Sd2r4   v <- rep(0,length=d^(2*r+4))  for(j in 0:(r+2))    v <- v+((-1)^j*OF(2*j)*choose(2*r+4, 2*j))*(K.pow(mu,2*r-2*j+4)%x%K.pow(vec(Sigma),j))    gamr<-(-1)^r*dmvnorm(mu,mean=rep(0,d),sigma=Sigma)*sum(w %*% v)  return(gamr)}################################################################################ Omega matrices (for exact MISE for normal mixtures)## Parameters # mus - means# Sigmas - variances# k - number of mixture components# a - subscript of Omega matrix# H - bandwidth matrix## Returns # Omega matrix###############################################################################omega <- function(mus, Sigmas, k, a, H, d, r, Sd2r){  ## the (i,j) element of Omega matrix is dmvnorm(0, mu_i - mu_j,  ## a*H + Sigma_i + Sigma_j)   if (k == 1)    omega.mat <- gamma.r(mu=rep(0,d),Sigma=a*H + 2*Sigmas, d=d, r=r, Sd2r=Sd2r)  ##dmvnorm(x=mus, mean=mus, sigma=a*H + 2*Sigmas)  else  {       if (is.matrix(mus)) d <- ncol(mus)    else d <- length(mus)    omega.mat <- matrix(0, nr=k, nc=k)    for (i in 1:k)    {      Sigmai <- Sigmas[((i-1)*d+1):(i*d),]      mui <- mus[i,]      for (j in 1:k)      {        Sigmaj <- Sigmas[((j-1)*d+1):(j*d),]        muj <- mus[j,]            omega.mat[i,j] <- gamma.r(mu=mui-muj, Sigma=a*H + Sigmai + Sigmaj, d=d, r=r, Sd2r=Sd2r) ## dmvnorm(x=mui, mean=muj, sigma=a*H + Sigmai + Sigmaj)      }    }  }    return(omega.mat)}omega.1d <- function(mus, sigmas, k, a, h, d=1, r, Sd2r){  ## the (i,j) element of Omega matrix is dmvnorm(0, mu_i - mu_j,  ## a*H + sigma_i + sigma_j)  H <- h^2  Sigmas <- sigmas^2    if (k == 1)    omega.mat <- gamma.r(mu=0, Sigma=as.matrix(a*H + 2*Sigmas), d=d, r=r, Sd2r=Sd2r)  ##dmvnorm(x=mus, mean=mus, sigma=a*H + 2*Sigmas)  else  {       omega.mat <- matrix(0, nr=k, nc=k)    for (i in 1:k)    {      Sigmai <- Sigmas[i]      mui <- mus[i]      for (j in 1:k)      {        Sigmaj <- Sigmas[j]        muj <- mus[j]            omega.mat[i,j] <- gamma.r(mu=mui-muj, Sigma=as.matrix(a*H + Sigmai + Sigmaj), d=d, r=r, Sd2r=Sd2r) ## dmvnorm(x=mui, mean=muj, sigma=a*H + Sigmai + Sigmaj)      }    }  }    return(omega.mat)}############################################################################### Exact MISE for normal mixtures## Parameters# mus - means# Sigmas - variances# Props - vector of proportions of each mixture component # H - bandwidth matrix# samp - sample size## Returns# Exact MISE for normal mixtures###############################################################################mise.mixt <- function(H, mus, Sigmas, props, samp, h, sigmas, deriv.order=0){  if (!(missing(h)))    return(mise.mixt.1d(h=h, mus=mus, sigmas=sigmas, props=props, samp=samp, deriv.order=deriv.order))     if (is.vector(mus)) d <- length(mus)  else d <- ncol(mus)  k <- length(props)  r <- deriv.order  Sd2r <- Sdr(d,2*r)  ## formula is found in Wand & Jones (1993) and Chacon, Duong & Wand (2008)  if (k == 1)   {    mise <- 2^(-r)*nu(r,H)/(samp * (4 * pi)^(d/2) * sqrt(det(H))) +         (1-1/samp)*omega(mus, Sigmas, 1, 2, H, d, r, Sd2r) -        2*omega(mus, Sigmas, 1, 1, H, d, r, Sd2r) +          omega(mus, Sigmas, 1, 0, H, d, r, Sd2r)  }  else  {    mise <- 2^(-r)*nu(r,H)/(samp * (4 * pi)^(d/2) * sqrt(det(H))) +      props %*% ((1-1/samp)*omega(mus, Sigmas, k, 2, H, d, r, Sd2r) -                  2*omega(mus, Sigmas, k, 1, H, d, r, Sd2r) +                  omega(mus, Sigmas, k, 0, H, d, r, Sd2r)) %*% props  }  return(drop(mise)) }mise.mixt.1d <- function(h, mus, sigmas, props, samp, deriv.order=0){    d <- 1  k <- length(props)  r <- deriv.order  Sd2r <- Sdr(d,2*r)  H <- as.matrix(h^2)    ## formula is found in Wand & Jones (1993) and Chacon, Duong & Wand (2008)  if (k == 1)   {    mise <- 2^(-r)*nu(r,H)/(samp * (4 * pi)^(d/2) * sqrt(det(H))) +         (1-1/samp)*omega.1d(mus, sigmas, 1, 2, h, d, r, Sd2r) -        2*omega.1d(mus, sigmas, 1, 1, h, d, r, Sd2r) +          omega.1d(mus, sigmas, 1, 0, h, d, r, Sd2r)  }  else  {    mise <- 2^(-r)*nu(r,H)/(samp * (4 * pi)^(d/2) * sqrt(det(H))) +      props %*% ((1-1/samp)*omega.1d(mus, sigmas, k, 2, h, d, r, Sd2r) -                  2*omega.1d(mus, sigmas, k, 1, h, d, r, Sd2r) +                  omega.1d(mus, sigmas, k, 0, h, d, r, Sd2r)) %*% props  }  return(drop(mise)) } ################################################################################ Exact AMISE for bivariate normal mixtures## Parameters# mus - means# Sigmas - variances# props - mixing proportions # H - bandwidth matrix# samp - sample size## Returns   # Exact AMISE for normal mixtures###############################################################################amise.mixt <- function(H, mus, Sigmas, props, samp, h, sigmas, deriv.order=0){  if (!(missing(h)))    return(amise.mixt.1d(h=h, mus=mus, sigmas=sigmas, props=props, samp=samp, deriv.order=deriv.order))  r <- deriv.order  if (is.vector(mus)) {d <- length(mus); mus <- t(matrix(mus))}  else d <- ncol(mus)  k <- length(props)   Sd2r4 <- Sdr(d,2*r+4)  ##w <- Sd2r4%*%(K.pow(vec(diag(d)),r)%x%K.pow(vec(H),2))  ##w <- as.vector(w)  if (k == 1)    omega.mat <- gamma.r2(mu=rep(0,d),Sigma=2*Sigmas, d=d, r=r, Sd2r4=Sd2r4, H=H)  else  {       omega.mat <- matrix(0, nr=k, nc=k)    for (i in 1:k)    {      Sigmai <- Sigmas[((i-1)*d+1):(i*d),]      mui <- mus[i,]      for (j in 1:k)      {        Sigmaj <- Sigmas[((j-1)*d+1):(j*d),]        muj <- mus[j,]            omega.mat[i,j] <- gamma.r2(mu=mui-muj, Sigma= Sigmai + Sigmaj, d=d, r=r, Sd2r4=Sd2r4, H=H)      }    }  }  if (k == 1) {    amise <- 2^(-r)*nu(r,H)/(samp * (4 * pi)^(d/2) * sqrt(det(H))) + omega.mat/4  }  else {    amise <- 2^(-r)*nu(r,H)/(samp * (4 * pi)^(d/2) * sqrt(det(H))) +      (props %*% omega.mat %*% props)/4  }   return(drop(amise))}amise.mixt.1d <- function(h, mus, sigmas, props, samp, deriv.order=0){  d <- 1  r <- deriv.order  k <- length(props)  H <- as.matrix(h^2)  Sd2r4 <- Sdr(d,2*r+4)  if (k == 1)    omega.mat <- gamma.r2(mu=rep(0,d),Sigma=2*sigmas^2, d=d, r=r, Sd2r4=Sd2r4, H=H)  else  {       omega.mat <- matrix(0, nr=k, nc=k)    for (i in 1:k)    {      Sigmai <- sigmas[i]^2      mui <- mus[i]      for (j in 1:k)      {        Sigmaj <- sigmas[j]^2        muj <- mus[j]            omega.mat[i,j] <- gamma.r2(mu=mui-muj, Sigma= Sigmai + Sigmaj, d=d, r=r, Sd2r4=Sd2r4, H=H)      }    }  }  if (k == 1) {    amise <- 2^(-r)*nu(r,H)/(samp * (4 * pi)^(d/2) * sqrt(det(H))) + omega.mat/4  }  else {    amise <- 2^(-r)*nu(r,H)/(samp * (4 * pi)^(d/2) * sqrt(det(H))) +      (props %*% omega.mat %*% props)/4  }   return(drop(amise))}amise.mixt.old <- function(H, mus, Sigmas, props, samp){   if (is.vector(mus)) {d <- length(mus); mus <- t(matrix(mus))}  else d <- ncol(mus)  k <- length(props)  Xi <- matrix(0, nr=k, nc=k)  for (i in 1:k)  {    Sigmai <- Sigmas[((i-1)*d+1) : (i*d),]    mui <- mus[i,]    for (j in 1:k)    {               Sigmaj <- Sigmas[((j-1)*d+1) : (j*d),]       muj <- mus[j,]       Aij <- chol2inv(chol(Sigmai + Sigmaj))       Bij <- Aij %*% (diag(d) - 2*(mui - muj) %*%  t(mui - muj) %*% Aij)       Cij <- Aij %*% (diag(d) - (mui - muj) %*%  t(mui - muj) %*% Aij)           Xi[i,j] <- dmvnorm.mixt(x=mui, mus=muj, Sigmas=Sigmai+Sigmaj, props=1) *                  (2*tr(H %*% Aij %*% H %*% Bij) + tr(H %*% Cij)^2)    }    }     amise <- 1/(samp *(4*pi)^(d/2)*sqrt(det(H)))+ 1/4*props %*% Xi %*% props  return(drop(amise))}################################################################################ Lambda matrices (for exact AMISE for normal mixtures)## Parameters # mus - means# Sigmas - variances