normal.r

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

  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),])
}