normal.r

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

###############################################################################
## 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 Sigma
dmvnorm.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 Sigma
dmvnorm.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)