prelim.r

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

  ##orep <- final.len <- prod(sapply(args, length))
  orep <- prod(sapply(args, length))
  
  for (i in 1:nargs) {
    x <- args[[i]]
    nx <- length(x)
    orep <- orep/nx
    cargs[[i]] <- rep(rep(x, rep(rep.fac, nx)), orep)
    rep.fac <- rep.fac * nx
  }
  do.call("cbind", cargs)
} 


###############################################################################
# For a list of matrices, this returns the dimensions of each matrix
###############################################################################

list.length <- function(x)
{
  ell <- length(x)
  len <- matrix(0, nr=ell, nc=2)
  for (i in 1:ell)
    len[i,] <- dim(x[[i]])

  return(len)
}

###############################################################################

permute.mat <- function(order)
{
    m <- as.integer(order)
    m <- m + 1
    eye <- diag(m)
    u1 <- eye[1:m,1]
    u2 <- eye[1:m,2:m]
    q1 <- kronecker(eye, u1)
    q2 <- kronecker(eye, u2)
    q <- matrix(c(q1, q2), nrow = nrow(q2), ncol = ncol(q1) + ncol(q2))
    if (!is.integer(q))
        storage.mode(q) <- "integer"
    q
}



###############################################################################
## Boolean functions
###############################################################################

is.even <- function(x)
{
  y <- x[x>0] %%2
  return(identical(y, rep(0, length(y))))
}

is.diagonal <- function(x)
{
  return(identical(diag(diag(x)),x))
}


####################################################################
## Functions for unconstrained pilot selectors, and (A)MISE-optimal
## selectors for normal mixtures
## Author: Jose E. Chacon
####################################################################


differences <- function(x, upper=TRUE)
{
  if (is.vector(x)) x <- t(as.matrix(x))
  n <- nrow(x)
  d <- ncol(x)
  
  difs <- matrix(ncol=d,nrow=n^2)
  for (j in 1:d)
  {    
    xj <- x[,j]
    difxj <- as.vector(xj%*%t(rep(1,n))-rep(1,n)%*%t(xj))
    ##The jth column of difs contains all the differences X_{ij}-X_{kj}
    difs[,j]<-difxj
  }
  
  if (upper)
  {
    ind.remove <- numeric()
    for (j in 1:(n-1))
      ind.remove <- c(ind.remove, (j*n+1):(j*n+j))
    
    return(difs[-ind.remove,])
  }
  else
    return(difs)
}


##### Odd factorial

OF<-function(m){factorial(m)/(2^(m/2)*factorial(m/2))} 


#####  Matrix square root

#matrix.sqrt <- function(A) {
#  d<-ncol(A)
#  if(d==1){A12<-sqrt(A)}
#  else{
#    xe <- eigen(A)
#    xe1 <- xe$values
#    if(all(xe1 >= 0)) {
#      xev1 <- diag(sqrt(xe1))
#    }
#    xval1 <- cbind(xe$vectors)
#    xval1i <- solve(xval1)
#    A12 <- xval1 %*% xev1 %*% xval1i}
#  return(A12)
#}


K.sum <- function(A,B)
{
  AB <- numeric()
  for (i in 1:nrow(A))
    for (j in 1:nrow(B))
      AB <- rbind(AB, A[i,] + B[j,])

  return(AB)
}


##### Commutation matrix of order m,n

K.mat<-function(m,n){
  K<-0
  for(i in 1:m){for(j in 1:n){
    H<-matrix(0,nrow=m,ncol=n)
    H[i,j]<-1
    K<-K+(H%x%t(H))
  }}
  return(K)        
}


##### Kronecker power of a matrix A

K.pow<-function(A,pow){
  if(floor(pow)!=pow)stop("pow must be an integer")
  Apow<-A
  if(pow==0){Apow<-1}
  if(pow>1){
    for(i in 2:pow) Apow<-Apow%x%A
  }
  return(Apow)
}
    
##### Row-wise Kronecker products and powers of matrices
    
mat.Kprod<-function(U,V){ #### Returns a matrix with rows U[i,]%x%V[i,]
  n1<-nrow(U)
  n2<-nrow(V)
  if(n1!=n2)stop("U and V must have the same number of vectors")
  p<-ncol(U)
  q<-ncol(V)
  P<-U[,1]*V
  if(p>1){for(j in 2:p){P<-cbind(P,U[,j]*V)}}
  return(P)
}

mat.K.pow<-function(A,pow){ #### Returns a matrix with the pow-th Kronecker power of A[i,] in the i-th row
  Apow<-A
  if(pow>1){
    for(i in 2:pow) Apow<-mat.Kprod(Apow,A)
  }
  return(Apow)
}

##### Vector of all r-th partial derivatives of the normal density at x=0, i.e., D^{\otimes r)\phi(0), for r=6,4

D6L0<-function(d,Sd6){
  r<-6
  DL0<-(-1)^(r/2)*(2*pi)^(-d/2)*OF(r)*(Sd6%*%K.pow(A=vec(diag(d)),pow=r/2))
  return(DL0)
}
    
D4L0<-function(d,Sd4){
  r<-4
  DL0<-(-1)^(r/2)*(2*pi)^(-d/2)*OF(r)*(Sd4%*%K.pow(A=vec(diag(d)),pow=r/2))
  return(DL0)
}

D2L0<-function(d,Sd2){
  r<-2
  DL0<-(-1)^(r/2)*(2*pi)^(-d/2)*OF(r)*(Sd2%*%K.pow(A=vec(diag(d)),pow=r/2))
  return(DL0)
}


T<-function(d,r){    #### Second version, recursive
  Id<-diag(d)
  Tmat<-Id
  Kdd<-K.mat(d,d)
  if(r>1){for(j in 2:r){
    Idj2Kdd<-diag(d^(j-2))%x%Kdd
    Tmat<-Idj2Kdd%*%((Tmat%x%Id)+Idj2Kdd)%*%Idj2Kdd
  }}
  return(Tmat)
}

## symmetriser matrix

Sdr<-function(d,r){
  if(r==0)S<-1
  else{
    S<-diag(d)
    if(r>=2){for(j in 2:r){S<-(S%x%diag(d))%*%T(d,j)/j}}}
  return(S)
}


##############################################################################
## Density derivative (psi) functional estimators 
##############################################################################


deriv.list <- function(d, r)
{
  derivt <- numeric()

  if (d==2)
  {
    for (j1 in r:0)
      for (j2 in r:0)
        if (sum(c(j1,j2))==r)
          derivt <- rbind(derivt, c(j1,j2))
  }
  if (d==3)
  {
    for (j1 in r:0)
      for (j2 in r:0)
        for (j3 in r:0)
          if (sum(c(j1,j2,j3))==r)
            derivt <- rbind(derivt, c(j1,j2,j3))
  }
  if (d==4)
  {
    for (j1 in r:0)
      for (j2 in r:0)
        for (j3 in r:0)
          for (j4 in r:0)
            if (sum(c(j1,j2,j3,j4))==r)
              derivt <- rbind(derivt, c(j1,j2,j3,j4))
  }

  if (d==5)
  {
    for (j1 in r:0)
      for (j2 in r:0)
        for (j3 in r:0)
          for (j4 in r:0)
            for (j5 in r:0)
              if (sum(c(j1,j2,j3,j4,j5))==r)
                derivt <- rbind(derivt, c(j1,j2,j3,j4,j5))
  }
  
  if (d==6)
  {
    for (j1 in r:0)
      for (j2 in r:0)
        for (j3 in r:0)
          for (j4 in r:0)
            for (j5 in r:0)
              for (j6 in r:0)
                if (sum(c(j1,j2,j3,j4,j5,j6))==r)
                  derivt <- rbind(derivt, c(j1,j2,j3,j4,j5,j6))
  }

  return(derivt)

}


RKfun <- function(r)
{
  if (r==0)
    val <- 1/(2*sqrt(pi))
  else if (r==1)
    val <- 1/(4*sqrt(pi))
  else if (r==2)
    val <- 3/(8*sqrt(pi)) 
  else if (r==3)
    val <- 15/(16*sqrt(pi))
  else if (r==4)
    val <- 105/(32*sqrt(pi))
  else if (r==5)
    val <- 945/(64*sqrt(pi))
  else if (r==6)
    val <- 10395/(128*sqrt(pi))
  else if (r==7)
    val <- 135135/(256*sqrt(pi))
  else if (r==8)
    val <- 2027025/(512*sqrt(pi))
  
  return(val)
}


########################################################################
### Identifying elements of Psi_4 matrix
########################################################################

Psi4.elem <- function(k, kprime, d)
{
  ind <- function(k, d)  
  {
    j <- 1
    dprime <- 1/2*d*(d+1)
    if (k > dprime) stop ("k is larger than d'")
    while (j < d & !(((j-1)*d -1/2*(j-2)*(j-1) < k) & (k <= j*d -1/2*j*(j-1))))
      j <- j+1
    i <- k - (j-1)*d + 1/2*j*(j-1)
    
    return(c(i,j))
  }

  ij <- ind(k, d)
  ei <- elem(ij[1],d)
  ej <- elem(ij[2],d)
  ipjp <- ind(kprime, d)
  eip <- elem(ipjp[1],d)
  ejp <- elem(ipjp[2],d)
  psi4.ind <- ei + eip + ej + ejp
  coeff <- (2 - (ij[1]==ij[2])) * (2 - (ipjp[1]==ipjp[2]))
  
  return (c(coeff, psi4.ind))
}


Psi4.list <- function(d)
{
  coeff <- vector()
  psifun <- vector()
  dprime <- 1/2*d*(d+1)
  for (k in 1:dprime)
    for (kprime in 1:dprime)
    {
      coeff <- c(coeff, Psi4.elem(k, kprime, d)[1])
      psifun <- rbind(psifun, Psi4.elem(k, kprime, d)[-1]) 
    }

  return(list(coeff=coeff, psi=psifun))  
}


default.gridsize <- function(d)
{
  if (d==1)
    gridsize <- 401
  else if (d==2)
    gridsize <- rep(151,d)
  else if (d==3)
    gridsize <- rep(51, d)
  else if (d==4)
    gridsize <- rep(21, d)

  return(gridsize)
}