mise.r

r软件 另一款可以计算核估计的软件包 需安装r软件

# k - number of mixture components
# r - derivative (r1, r2)
#
# Returns 
# Lambda matrix
###############################################################################

lambda <- function(mus, Sigmas, k, r)
{
  ## the (i,j) element of Lambda matrix is d^r/ dx^r  dmvnorm(0, mu_i - mu_j,
  ## a*H + Sigma_i + Sigma_j)
    
  if (is.vector(mus)) d <- length(mus)
  else d <- ncol(mus)
  
  if (k == 1) 
    lambda.mat <- dmvnorm.deriv.2d(r=r, x=rep(0, length(mus)), Sigma=2*Sigmas)
  else
  {   
    if (is.matrix(mus)) d <- ncol(mus)
    else d <- length(mus)
    lambda.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,]    
        lambda.mat[i,j] <- dmvnorm.deriv.2d(r=r, x=mui-muj,Sigma=Sigmai+Sigmaj)
      }
    }
  }
  
  return(lambda.mat)
}


amise.mixt.2d <- function(H, mus, Sigmas, props, samp)
{  
  d <- ncol(Sigmas)
  k <- length(props)
  
  ##h1 <- sqrt(H[1,1])
  ##h2 <- sqrt(H[2,2])
  ##h12 <- H[1,2]

  ## formula is found in Wand & Jones (1993)
  if (k == 1) 
  {
   amise <- 1/(samp * (4*pi)^(d/2) * sqrt(det(H))) +
       1/4 *(lambda(mus, Sigmas, k, r=c(4,0))*H[1,1]^2 +
           4*lambda(mus, Sigmas, k, r=c(3,1))*H[1,1]*H[1,2] +  
           2*lambda(mus, Sigmas, k, r=c(2,2))*(H[1,1]*H[2,2] + 2*H[1,2]^2) + 
           4*lambda(mus, Sigmas, k, r=c(1,3))*H[2,2]*H[1,2]+    
             lambda(mus, Sigmas, k, r=c(0,4))*H[2,2]^2) 
  }
  else
  {
    amise <- 1/(samp * (4*pi)^(d/2) * sqrt(det(H))) +
      1/4 * props %*% 
          (  lambda(mus, Sigmas, k, r=c(4,0))*H[1,1]^2 +
           4*lambda(mus, Sigmas, k, r=c(3,1))*H[1,1]*H[1,2] +  
           2*lambda(mus, Sigmas, k, r=c(2,2))*(H[1,1]*H[2,2] + 2*H[1,2]^2) + 
           4*lambda(mus, Sigmas, k, r=c(1,3))*H[2,2]*H[1,2]+    
             lambda(mus, Sigmas, k, r=c(0,4))*H[2,2]^2) %*% props
  }
  
  return(drop(amise)) 
}


###############################################################################
# Finds the bandwidth matrix that minimises the MISE for normal mixtures
#
# Parameters
# mus - means
# Sigmas - variances
# props - vector of proportions of each mixture component 
# Hstart - initial bandwidth matrix
# samp - sample size
# full - 1 minimise over full bandwidth matrices
#      - 0 minimise over diagonal bandwidth matrices
# 
# Returns
# H_MISE
###############################################################################

hmise.mixt <- function(mus, sigmas, props, samp, hstart, deriv.order=0)
{
  r <- deriv.order
  d <- 1
  
  if (missing(hstart))
  {
    x <- rnorm.mixt(n=1000, mus=mus, sigmas=sigmas)  
    hstart <- sqrt((4/(samp*(d+2*r+2)))^(2/(d+2*r+4)) * var(x))
  }
  mise.mixt.temp <- function(h)
  {  
    return(mise.mixt.1d(h=h, mus=mus, sigmas=sigmas, props=props, samp=samp, deriv.order=deriv.order))
  }

  result <- optimize(f=mise.mixt.temp, interval=c(0, 10*hstart))
  hmise <- result$minimum
  return(hmise)
}

Hmise.mixt <- function(mus, Sigmas, props, samp, Hstart, deriv.order=0)
{
  r <- deriv.order
  if (is.vector(mus)) d <- length(mus)
  else d <- ncol(mus) 
 
  ## use normal reference estimate as initial condition
  if (missing(Hstart))
  {
    x <- rmvnorm.mixt(10000, mus, Sigmas, props)
    Hstart <- matrix.sqrt((4/(samp*(d+2*r+2)))^(2/(d+2*r+4)) * var(x))
  }
  
  Hstart <- vech(Hstart)

  # input vech(H) into mise.mixt.temp because optim can only optimise
  # over vectors and not matrices
  mise.mixt.temp <- function(vechH)
  {  
    H <- invvech(vechH) %*% invvech(vechH)
    ## using H <- invvech(vechH) %*% invvech(vechH) ensures that H
    ## is positive definite
    
    return(mise.mixt(H=H, mus=mus, Sigmas=Sigmas, props=props, samp=samp, deriv.order=deriv.order))
  }

  result <- optim(Hstart, mise.mixt.temp, method="BFGS")
  Hmise <- invvech(result$par) %*% invvech(result$par) 
  
  return(Hmise)
}   

Hmise.mixt.diag <- function(mus, Sigmas, props, samp, Hstart, deriv.order=0)
{   
  if (is.vector(mus)) d <- length(mus)
  else d <- ncol(mus) 

  if (missing(Hstart))
  {
    x <- rmvnorm.mixt(10000, mus, Sigmas, props)
    Hstart <- (4/(samp*(d + 2)))^(2/(d + 4)) * var(x)
  }  
  Hstart <- diag(matrix.sqrt(Hstart))

  mise.mixt.temp <- function(diagH)
  {  
    H <- diag(diagH) %*% diag(diagH)
    return(mise.mixt(H=H, mus=mus, Sigmas=Sigmas, props=props, samp=samp, deriv.order=deriv.order))
  }

  result <- optim(Hstart, mise.mixt.temp, method = "BFGS")
  Hmise <- diag(result$par) %*% diag(result$par) 
  
  return(Hmise)
}   



###############################################################################
## Finds bandwidth matrix that minimises the AMISE for normal mixtures - 2-dim
##
## Parameters
## mus - means
## Sigmas - variances
## props - vector of proportions of each mixture component 
## Hstart - initial bandwidth matrix
## samp - sample size
## 
## Returns
## Bandwidth matrix that minimises AMISE
###############################################################################

hamise.mixt <- function(mus, sigmas, props, samp, hstart, deriv.order=0)
{
  r <- deriv.order
  d <- 1
 
  if (missing(hstart))
  {
    x <- rnorm.mixt(n=1000, mus=mus, sigmas=sigmas)  
    hstart <- sqrt((4/(samp*(d+2*r+2)))^(2/(d+2*r+4)) * var(x))
  }
  amise.mixt.temp <- function(h)
  {  
    return(amise.mixt.1d(h=h, mus=mus, sigmas=sigmas, props=props, samp=samp, deriv.order=deriv.order))
  }

  result <- optimize(f=amise.mixt.temp, interval=c(0, 10*hstart))
  hamise <- result$minimum
  
  return(hamise)
}


Hamise.mixt <- function(mus, Sigmas, props, samp, Hstart, deriv.order=0)
{
  r <- deriv.order
  if (is.vector(mus)) d <- length(mus)
  else d <- ncol(mus) 

  ## use explicit formula for single normal
  if (length(props)==1)
  {
    Hamise <- (4/ (samp*(d+2*r+2)))^(2/(d+2*r+4)) * Sigmas
  }
  else
  {  
    ## use normal reference estimate as initial condition
    if (missing(Hstart)) 
    {
      x <- rmvnorm.mixt(10000, mus, Sigmas, props)
      Hstart <- matrix.sqrt((4/ (samp*(d+2*r+2)))^(2/(d+2*r+4)) * var(x))
    }
    
    ## input vech(H) into mise.mixt.temp because optim can only optimise
    ## over vectors and not matrices    
    Hstart <- vech(Hstart)
    amise.mixt.temp <- function(vechH)
    {
      H <- invvech(vechH) %*% invvech(vechH)
      ## ensures that H is positive definite
      
      return(amise.mixt(H=H, mus=mus, Sigmas=Sigmas, props=props, samp=samp, deriv.order=deriv.order))
    }
    
    result <- optim(Hstart, amise.mixt.temp, method="BFGS")
    Hamise <- invvech(result$par) %*% invvech(result$par) 
  }
  
  return(Hamise)
}   
  

###############################################################################
# ISE for normal mixtures (fixed KDE)
# 
# Parameters
# x - data values
# H - bandwidth matrix
# 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 - mixing proportions
# lower - vector of lower end points of rectangle
# upper - vector of upper end points of rectangle
# gridsize - vector of number of grid points
# stepsize - vector of step sizes
# Returns
# ISE 
###############################################################################

ise.mixt <- function(x, H, mus, Sigmas, props, h, sigmas, deriv.order=0)
{
  if (!(missing(h)))
    return(ise.mixt.1d(x=x, h=h, mus=mus, sigmas=sigmas, props=props, deriv.order=deriv.order))
  
  ##if (is.list(x))
  ##  return (ise.mixt.pc(x, H, mus, Sigmas, props, lower, upper, gridsize,
  ##                      stepsize))
  if (is.vector(x)) x <- matrix(x,nr=1)
  if (is.vector(mus)) mus <- matrix(mus, nr=length(props))

  d <- ncol(x)
  n <- nrow(x)
  M <- length(props)
  ise1 <- 0
  ise2 <- 0
  ise3 <- 0

  ## formula is found in thesis  
  if (d==2)
    ise1 <- dmvnorm.2d.sum(x=x, Sigma=2*H, inc=1)
  else if (d==3)
    ise1 <- dmvnorm.3d.sum(x=x, Sigma=2*H, inc=1)
  else if (d==4)
    ise1 <- dmvnorm.4d.sum(x=x, Sigma=2*H, inc=1)
  else if (d==5)
    ise1 <- dmvnorm.5d.sum(x=x, Sigma=2*H, inc=1)
  else if (d==6)
    ise1 <- dmvnorm.6d.sum(x=x, Sigma=2*H, inc=1)
  
  for (j in 1:M)
  {
    Sigmaj <- Sigmas[((j-1)*d+1) : (j*d),]
    ise2 <- ise2 + sum(props[j]*dmvnorm(x=x, mean=mus[j,], sigma=H + Sigmaj))
    
    for (i in 1:M)
    {
      Sigmai <- Sigmas[((i-1)*d+1) : (i*d),]
      ise3 <- ise3 + sum(props[i] * props[j] *
                         dmvnorm(x=mus[i,], mean=mus[j,], sigma=Sigmai+Sigmaj))
    }
  }  

  return (ise1/n^2 - 2*ise2/n + ise3)
}


ise.mixt.1d <- function(x, h, mus, sigmas, props, deriv.order=0)
{  
  d <- 1
  n <- length(x)
  M <- length(props)
  ise1 <- 0
  ise2 <- 0
  ise3 <- 0
  
  ise1 <- dnorm.sum(x=x, sigma=sqrt(2)*h, inc=1)
  
  for (j in 1:M)
  {
    sigmaj <- sigmas[j]
    ise2 <- ise2 + sum(props[j]*dnorm(x=x, mean=mus[j], sd=sqrt(h^2 + sigmaj^2)))
    
    for (i in 1:M)
    {
      sigmai <- sigmas[i]
      ise3 <- ise3 + sum(props[i]*props[j]*dnorm(x=mus[i], mean=mus[j], sd=sqrt(sigmai^2+sigmaj^2)))
    }
  }  

  return (ise1/n^2 - 2*ise2/n + ise3)
}