prelim.r

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

#############################################################################
## Basic important vectors and matrices
#############################################################################
  
###############################################################################
# Vec operator 
#
# Takes square matrix x and returns its elements stacked column wise in a
# vector
###############################################################################

vec <- function(x, byrow=FALSE)
{
  if (is.vector(x)) return (x)
  
  if (byrow) x <- t(x)
  d <- ncol(x)
  vecx <- vector()
  for (j in 1:d)
    vecx <- c(vecx, x[,j])
  
  return(vecx)           
}

###############################################################################
# Vech operator 
# 
# Takes matrix x and returns its elements that in the lower triangular half, 
# stacked column wise in a vector
###############################################################################

vech <- function(x)
{
  if (is.vector(x))
  {
    if (length(x)==1)
      return (x)
    else
      stop("vech undefined for vectors")
  }
  else if (is.matrix(x))
  {  
    d <- ncol(x)
    if (d!=nrow(x)) ##if (!isSymmetric(x))
      stop("vech only defined for square matrices")
    
    vechx <- vector()
    for (j in 1:d)
      vechx <- c(vechx, x[j:d,j])
    return(vechx)
  }
  
}

# Analogue for stacked matrix
vech.cat <- function(x)
{
  d <- ncol(x)
  num <- nrow(x)/ncol(x)
  vechx <- vector()
  
  for (j in 1:num)
    vechx <- c(vechx, vech(x[((j-1)*d+1) : (j*d),]))
                   
  return(vechx)
}

  
###############################################################################
# Inverse vec operator 
#
# Takes vector x and stacks its elements into columns of a matrix. 
###############################################################################
    
invvec <- function(x, ncol, nrow, byrow=FALSE)
{
  if (length(x)==1)
    return(x)
  
  d <- sqrt(length(x))
  if (missing(ncol) | missing(nrow))
  {
    ncol <- d; nrow <- d
    if (round(d) != d)
      stop("Need to specify nrow and ncol for non-square matrices")
  }
  
  invvecx <- matrix(0, nrow = nrow, ncol = ncol)
  if (byrow)
    for (j in 1:nrow)
      invvecx[j,] <- x[c(1:ncol) + (j-1)*ncol]
  else
    for (j in 1:ncol)
      invvecx[,j] <- x[c(1:nrow) + (j-1)*nrow]
  
  return(invvecx)
}

###############################################################################
# Inverse vech operator 
#
# Takes vector x and stacks its elements into a symmetric matrix 
###############################################################################

invvech <- function(x)
{
  if (length(x)==1)
    return(x)
  
  d <- (-1 + sqrt(8*length(x) + 1))/2
  if (round(d) != d)
    stop("Number of elements in x will not form a square matrix.")
  invvechx <- matrix(0, nr = d, nc = d)

  for (j in 1:d)
    invvechx[j:d,j] <- x[1:(d-j+1)+ (j-1)*(d - 1/2*(j-2))]
  
  invvechx <- invvechx + t(invvechx) - diag(diag(invvechx))
  
  return(invvechx)
}

##### Trace of matrix

tr <- function(A)
{
  count <- 0
  if (is.vector(A)) return (A[1])
  if (nrow(A)!=ncol(A))
    stop('Not square matrix')

  else 
     for (i in 1:nrow(A))
       count <- count + A[i,i]
  return(count)
}


###############################################################################
# Elementary vector 
# 
# Creates a vector of length d and with 1 at i-th component and 0 elsewhere 
###############################################################################
    
elem <- function(i, d)
{
  elem.vec <- rep(0, d)
  elem.vec[i] <- 1
  
  return(elem.vec)
}      

###############################################################################
# Matrix square root - taken from Stephen Lake 
# http://www5.biostat.wustl.edu/s-news/s-news-archive/200109/msg00067.html
###############################################################################

matrix.sqrt <- function(A)
{
  if (length(A)==1)
    return(sqrt(A))
  sva <- svd(A)
  if (min(sva$d)>=0)
    Asqrt <- sva$u %*% diag(sqrt(sva$d)) %*% t(sva$v)
  else
    stop("Matrix square root is not defined")
  return(Asqrt)
}

###############################################################################
# Duplication matrix
# Taken from Felipe Osorio http://www.ime.usp.br/~osorio/files/dupl.q
###############################################################################

dupl <- function(order, ret.q = FALSE)
{
    # call
    cl <- match.call()
    time1 <- proc.time()
    if (!is.integer(order))
        order <- as.integer(order)
    n <- order - 1
    
    # initial duplication matrix
    d1 <- matrix(0, nrow = 1, ncol = 1)
    d1[1,1] <- 1
    if (!is.integer(d1))
        storage.mode(d1) <- "integer"
    
    # recursive formula
    if (n > 0){
    	for (k in 1:n){
    	    drow <- 2*k + 1 + nrow(d1)
    	    dcol <- k + 1 + ncol(d1)
    	    d2 <- matrix(0, nrow = drow, ncol=dcol)
    	    storage.mode(d2) <- "integer"
    	    d2[1,1] <- 1
    	    d2[2:(k+1),2:(k+1)] <- diag(k)
    	    d2[(k+2):(2*k+1),2:(k+1)] <- diag(k)
    	    d2[(2*k+2):drow,(k+2):dcol] <- d1
    	    # permutation matrix
    	    q <- permute.mat(k)
    	    # new duplication matrix
    	    d2 <- q %*% d2
    	    storage.mode(d2) <- "integer"
    	    d1 <- d2
    	}
    }
    else {
    	d2 <- q <- d1
    }
    
    # results
    obj <- list(call=cl, order=order, d=d2)
    if (ret.q)
        obj$q <- q
    obj$time <- proc.time() - time1
    obj
}

invdupl <- function(order, ret.q = FALSE)
{
    # call
    cl <- match.call()
    time1 <- proc.time()
    if (!is.integer(order))
        order <- as.integer(order)
    n <- order - 1

    # initial inverse of duplication matrix
    h1 <- matrix(0, nrow = 1, ncol = 1)
    h1[1,1] <- 1

    # recursive formula
    if (n > 0){
    	for (k in 1:n){
    	    hrow <- k + 1 + nrow(h1)
    	    hcol <- 2*k + 1 + ncol(h1)
    	    h2 <- matrix(0, nrow = hrow, ncol=hcol)
    	    h2[1,1] <- 1
    	    h2[2:(k+1),2:(k+1)] <- .5*diag(k)
    	    h2[2:(k+1),(k+2):(2*k+1)] <- .5*diag(k)
    	    h2[(k+2):hrow,(2*k+2):hcol] <- h1
    	    # permutation matrix
    	    q <- permute.mat(k)
    	    # new inverse of duplication matrix
    	    h2 <- h2 %*% t(q)
    	    h1 <- h2
    	}
    }
    else {
    	h2 <- q <- h1
    }
    
    # results
    obj <- list(call=cl, order=order, h=h2)
    if (ret.q)
        obj$q <- q
    obj$time <- proc.time() - time1
    obj
}

###############################################################################
# Pre-sphering
# Parameters
# x - data points
#
# Returns
# Pre-sphered x values
###############################################################################

pre.sphere <- function(x, mean.centred=FALSE)
{
  S <- var(x)
  Sinv12 <- matrix.sqrt(chol2inv(chol(S)))

  if (mean.centred)
  {
    xmean <- apply(x,2,mean)
    for (i in 1:ncol(x))
      x[,i] <- x[,i] - xmean[i]
  }
  x.sphered <- matrix(0, nc=ncol(x), nr=nrow(x))
  for (i in 1:nrow(x))
    x.sphered[i,] <- Sinv12 %*% x[i,]    

  return (x.sphered)
}

pre.sphere.pc <- function(x.pc)
{
  g <- length(x.pc$nclust)
  d <- ncol(x.pc$x)
  x.pc1 <- x.pc
  x.pc1$x <- matrix(0, nc=d, nr=nrow(x.pc$x))
  for (j in 1:g)
  {
    xj <- x.pc$x[x.pc$ind==j,]
    x.pc1$x[x.pc$ind==j,] <- pre.sphere(xj) 
  }
 
  return (x.pc1)          
}  

###############################################################################
# Pre-scaling
# Parameters
# x - data points
#
# Returns
# Pre-scaled x values
###############################################################################

pre.scale <- function(x, mean.centred=FALSE)
{
  x.scaled <- numeric()
  x.sd <- apply(x, 2, sd)
  d <- ncol(x)

  for (i in 1:d)
    if (mean.centred)
      x.scaled <- cbind(x.scaled, (x[,i] - mean(x[,i]))/x.sd[i])
    else
      x.scaled <- cbind(x.scaled, x[,i]/x.sd[i])                  

  return (x.scaled)
}


###############################################################################
# Finds row index matrix
# Parameters
# x - data points
#
# Returns
# i  - if r==mat[i,]
# NA - otherwise
###############################################################################
which.mat <- function(r, mat)
{
  ind <- numeric()
  
  for (i in 1:nrow(mat))
    if (identical(r, mat[i,])) ind <- c(ind,i)

  return(ind)  
}



###############################################################################
# Permute a list of values
#
# Same function as EXPAND.GRID (base package), modified to take 
# list as an argument and returns a matrix 
###############################################################################

permute <- function (args) 
{
  nargs <- length(args)
  if (!nargs) 
    return(as.data.frame(list()))
  if (nargs == 1 && is.list(a1 <- args[[1]])) 
    nargs <- length(args <- a1)
  if (nargs <= 1) 
    return(as.data.frame(if (nargs == 0 || is.null(args[[1]])) list() else args, 
                         optional = TRUE))
  cargs <- args
  rep.fac <- 1