多次元ベクトルデータの超球面表現と密度推定2

• Rのkde2d関数の多次元版…スケーリングとか、bandwidth.nrd関数からの出力の補正とかの点でちょっと怪しい…

kdeNd<-function (x,h, n = 25) 
{
    max<-apply(x,2,max)
    min<-apply(x,2,min)
    nx <- ncol(x)#No. axis
    ns <- nrow(x)#No. samples
    g<-mapply(seq.int,min,max,length.out=n)

    if (missing(h)) 
        h <- apply(x,2,bandwidth.nrd)
    #h <- h/4
    h <-h/2^nx
    
    a<-list(NULL)
    
    for(i in 1:nx){
     a[[i]]<-outer(g[,i],x[,i],"-")/h[i]
    }
    ret<-list(NULL)
    counter<-0
    kurai<-rep(1,nx)
    N<-n^nx
    for(i in 1:N){


     for(j in 1:(length(kurai)-1)){
      if(kurai[j]==n+1){
       kurai[j]<-1
       kurai[j+1]<-kurai[j+1]+1
      }
     }
     ret[[counter<-counter+1]]<-c(rep(0,nx+1))
     #ret[[counter]][1]<-1/ns
     tmp<-c(rep(1,ns))
     tmp2<-0
     for(j in 1:nx){
      
      for(k in 1:ns){
       tmp[k]<-tmp[k]*dnorm(a[[j]][kurai[j],k])

      }
     }
     for(j in 1:ns){
      tmp2<-tmp2+tmp[j]
     }
     for(j in 1:nx){
     
      tmp2<-tmp2/h[j]
      
      ret[[counter]][j+1]<-g[kurai[j],j]
     }
     ret[[counter]][1]<-tmp2/(ns^(nx-1))
     kurai[1]<-kurai[1]+1
    }

    #z <- matrix(dnorm(ax), n, nx) %*% t(matrix(dnorm(ay), n, nx))/(nx * h[1] * h[2])
    return(ret)
}

• ちなみにkde2dは

function (x, y, h, n = 25, lims = c(range(x), range(y))) 
{
    nx <- length(x)
    if (length(y) != nx) 
        stop("data vectors must be the same length")
    if (any(!is.finite(x)) || any(!is.finite(y))) 
        stop("missing or infinite values in the data are not allowed")
    if (any(!is.finite(lims))) 
        stop("only finite values are allowed in 'lims'")
    gx <- seq.int(lims[1], lims[2], length.out = n)
    gy <- seq.int(lims[3], lims[4], length.out = n)
    if (missing(h)) 
        h <- c(bandwidth.nrd(x), bandwidth.nrd(y))
    h <- h/4
    ax <- outer(gx, x, "-")/h[1]
    ay <- outer(gy, y, "-")/h[2]
    z <- matrix(dnorm(ax), n, nx) %*% t(matrix(dnorm(ay), n, 
        nx))/(nx * h[1] * h[2])
    return(list(x = gx, y = gy, z = z))
}
<environment: namespace:MASS>