カーネル推定のオーダー

  • こちらでGaussian Sequence Modelをなぞっている。その一環
  • カーネル関数にはオーダーというものがあって、畳みこんでスムージングするときのバイアス制御と関係するという話があり、それを用いて話が進んで行くのだが、ぼんやりしているので、確認しておく
  • 参考はこちら
  • まず、カーネル関数は、台全体での積分が1である。確率密度関数とこの点は同様だが、確率密度関数が非負であるのに対して、カーネル関数の値は負でもよい
    • \int_{-\infty}^{\infty} k(v) dv = 1
  • カーネル関数の0周りのp次モーメント\mu_p=\int_{-\infty}^{\infty} v^p k(v) dvによってカーネル関数を特徴分類する
  • 上記のように\mu_0 = 1
  • (後述のように・後述しないかもしれないけれど)カーネル関数のモーメントは、カーネル推定において、バイアスと関連しているので、0であるモーメントが多いと便利である
  • そのため、\mu_1 = 1であるようなカーネルは大事。これは、「平均が0」ということ。逆に言うと、これはバイアス制御上大事なので、0周りのp次モーメントを考える代わりに\mu_1周りのモーメントを考えることが多い。それはカーネル関数の形を変えずにシフトするだけ。したがって、カーネル関数の形に本質的に影響しない
  • それ以外でp次モーメントに0のものを持ち出すには、0周り(\mu_1=0周り)で対称な関数を考えることである。対称な関数では奇数次モーメントがすべて0になる。したがって、対称関数は都合がよい
  • というわけで対称なカーネル関数を考える
  • 対称なカーネル関数の奇数次モーメントは0であるから、さらなる分類をするとなると、偶数次モーメントが0か否かに注目することになる。
  • q-thオーダーモーメント関数と言うとき、q未満のすべてのモーメントが0であるような関数であると言うことにする
  • ちなみに、非負の値しかとらないカーネル関数では、\mu_2 > 0であるから、必ず2nd order カーネルとなる
  • 2nd order カーネルの例
    • Uniform k_0(v) = \frac{1}{2} \text{when } |v| \le 1
    • Epanchnikov k_1(v) = \frac{3}{4}(1-u^2) \text{when } |v| \le 1
    • Boweight k_2(v) = \frac{15}{16}(1-u^2)^2 \text{when } |v| \le 1
    • Triweight k_3(v) = \frac{35}{32}(1-u^2)^3 \text{when } |v| \le 1
    • Gaussian k_{gauss}(v) = \frac{1}{\sqrt{2 \pi}}exp^{-\frac{v^2}{2}}

  • じゃあ4-th orderカーネルは…
    • 2nd orderカーネルの周辺部に負の領域を抱き込ませた形
    • Epanchnikov k_{4,1}(v) = \frac{15}{8})(1-\frac{7}{3}v^2)\times k_1(v)
    • Biweight k_{4,2}(v) = \frac{7}{4}(1-3v^2)\times k_2(v)
    • Triweight k_{4,3}(v) = \frac{27}{16}(1-\frac{11}{3}u^2)\times k_3(v)
    • Gaussian k_{4,gauss}(v) = \frac{1}{2}(3-u^2)\times k_{gauss}(v)

  • さらに6-th orderカーネルは…
    • Epanechnikov k_{6,1}(v) = \frac{175}{64}(1-6v^2+\frac{33}{5}v^4)k_1(v)
    • Biweight k_{6,2}(v) = \frac{315}{128}(1-\frac{22}{3}v^2+\frac{143}{15}v^4)k_2(v)
    • Triweight k_{6,3}(v) = \frac{297}{128}(1-\frac{26}{3}v^2+13v^2)k_3(v)
    • Gaussian k_{6,gauss}(v) = \frac{1}{8}(15-10v^2+v^4)k_{gauss}(v)

  • Gaussian-basedカーネルの2nd, 4th, 6th オーダーカーネルを重ねてプロットしてみる
    • 実際、Gaussian-based カーネルの2r-th order関数は、2nd order gaussian kernel (k_{gauss})の導関数を用いて以下のような形をして表現される(参考こちら)
      • k_{2r,gauss}=\frac{(-1)^r*k_{gauss}^{(2r-1)}(v)}{2^{r-1}(r-1)!x}

k0 <- function(u){ # Uniform
	ret <- rep(0,length(u))
	a <- which(abs(u)<=1)
	ret[a] <- 1/2
	ret
}
k1 <- function(u){ # Epanchnikov
	ret <- rep(0,length(u))
	a <- which(abs(u)<=1)
	ret[a] <- 3/4 * (1-u[a]^2)
	ret
}
k2 <- function(u){ # Biweight
	ret <- rep(0,length(u))
	a <- which(abs(u)<=1)
	ret[a] <- 15/16 * (1-u[a]^2)^2
	ret
}
k3 <- function(u){ # Triweight
	ret <- rep(0,length(u))
	a <- which(abs(u)<=1)
	ret[a] <- 35/32 * (1-u[a]^2)^3
	ret
}
k.gauss <- function(u){ # Gaussian
	ret <- 1/sqrt(2*pi) * exp(-u^2/2)
	ret
}

u <- seq(from=-5,to=5,length=1000)
X.order2 <- matrix(0,length(u),5)
X.order2[,1] <- k0(u)
X.order2[,2] <- k1(u)
X.order2[,3] <- k2(u)
X.order2[,4] <- k3(u)
X.order2[,5] <- k.gauss(u)

matplot(u,X.order2,type="l")

k4.1 <- function(u){
	15/8*(1-7/3*u^2)*k1(u)
}
k4.2 <- function(u){
	7/4*(1-3*u^2)*k2(u)
}
k4.3 <- function(u){
	27/16*(1-11/3*u^2)*k3(u)
}
k4.gauss <- function(u){
	1/2*(3-u^2)*k.gauss(u)
}

X.order4 <- matrix(0,length(u),4)
X.order4[,1] <- k4.1(u)
X.order4[,2] <- k4.2(u)
X.order4[,3] <- k4.3(u)
X.order4[,4] <- k4.gauss(u)
matplot(u,X.order4,type="l")

k6.1 <- function(u){
	175/64*(1-6*u^2+33/5*u^4)*k1(u)
}
k6.2 <- function(u){
	315/128*(1-22/3*u^2+143/15*u^4)*k2(u)
}
k6.3 <- function(u){
	297/128*(1-26/3*u^2+13*u^4)*k3(u)
}
k6.gauss <- function(u){
	1/8*(15-10*u^2+u^4)*k.gauss(u)
}

X.order6 <- matrix(0,length(u),4)
X.order6[,1] <- k6.1(u)
X.order6[,2] <- k6.2(u)
X.order6[,3] <- k6.3(u)
X.order6[,4] <- k6.gauss(u)
matplot(u,X.order6,type="l")


matplot(u,cbind(X.order2[,5],X.order4[,4],X.order6[,4]),type="l")

カーネル推定と特性関数とモーメントとフーリエ変換

  • こちらでGaussian Sequence Modelをなぞっている。その一環
  • カーネル推定は、積分して1になる関数(カーネル関数)を重みづけ関数として、「元の関数(や観測データ)」を畳みこんで推定値(推定関数)を作成する方法
  • カーネル関数は(負の値を持つ場合もあるが)、確率密度関数のようなものである
  • 確率密度関数には、積率母関数(moment generating function)がある(場合がある)
    • M_X(t) := E(e^{tX})と書かれるが、「確率変数Xにこのように定義される関数」のこと
    • 積率母関数のt=0の周囲のn階導関数は、E(X^n) = M^{(n)}_X(0) = \frac{d^n M_X}{dt^n}(0)のように、Xのn次モーメントの期待値が積率母関数Mのn階微分の関数のt=0での式・値になる
  • この積率母関数は、確率密度関数そのものと同様に、確率変数をきちんと定義づける関数であるが、そのような性質を持つ関数に特性関数がある
  • 確率密度関数にある特性関数はカーネル関数にも定義できる(多分)。そうすると、モーメントの期待値が特性関数のt=0周りの導関数として取り扱える
  • カーネル関数がq-th order カーネルというとき、p=1,2,...,q-1において\int_{-\infty}^{\infty} v^p K(v) dv=0となる性質があり、このようなpが大きいとカーネル関数を用いた畳みこみ計算は楽(収束が速い?)になる。
    • この\int_{-\infty}^{\infty} v^p K(v) dv=0がモーメントのようなもので、したがって、カーネル関数のモーメントがどうなっているのか、q-th order カーネル関数である、というときのqがわかるとよい。
  • このqを知るのに、カーネル関数をフーリエ変換してやると、t=0周りのn階微分値を取り出すことは容易であって、カーネル関数の畳みこみにおける挙動の重要な情報が取り出せる
  • これが、カーネル推定と特性関数とモーメントとフーリエ変換フーリエ展開との関係
  • 実際には何をminimax推定するかというと、カーネル関数はK_h(u)=\frac{1}{h}K(\frac{u}{h})のように「hという幅が可変」であり、この値として「最適化」をしたいものであるから、それがminimax推定できる枠組みが得られればよい
  • カーネル関数で畳みこんで、それについて台全体に渡って「ずれ」の積分をすることが「リスク関数値」の算出になる
  • その計算をするにあたり、テイラー展開をかませると、台の着目値周りにおけるモーメント・導関数が効いてくる。どのくらいのオーダーまで考慮するかで、テイラー展開の精度が影響を受けるが、q-th order関数を用いて、q次まで見てやることにすれば、q-1次までのモーメントが0だから、台全体での積分によって消える(vanishing moments)。残るはq次だけ(それより高い次数は近似という意味合いで省略してある)。結局、ある程度の次数でのテイラー展開という近似をするときに、1,2,...,q次の項のすべてを計算せずに第q次だけの計算でよくなって便利、というのが、q-th orderカーネルの役どころ(らしい)
  • このようにして計算が楽になったリスクには、カーネル関数の幅を決めるパラメタh(K_h(u)=\frac{1}{h}K(\frac{u}{h}))が登場し、このhを小さくすると当てはまりは良くなるが、分散が大きくなる、というminimaxの構図に持ち込まれる

量的データで決断

  • 帰結が2種類あるときのオプション選択のことを考えた(こちらとか)
  • 帰結が3以上種類あるときのオプション選択のことを考えた(こちらとか)
  • 自然な流れとして、量的帰結についてのオプション選択を考えたい
  • 帰結の分布を観測データからオプション別に思い描いて、その分布に対して評価を下したい
  • まずは、帰結分布を描くことにする
  • カテゴリカルな帰結の場合は、共役事前分布であるベータ分布・ディリクレ分布を用いた
  • 量的な場合には、ガウシアンでの密度推定がよいだろう
    • 密度推定するには、平滑化をコントロールする必要があるが、それについては、リスク関数の最小化など、よく使うもの(Rのdensity()関数でのデフォルト)を決め打ちにしてみることにする

  • 書き散らしソース
true.prob <- list()
true.prob[[1]] <- matrix(c(1,1,0.1,10,2,0.9),byrow=TRUE,ncol=3)
true.prob[[2]] <- matrix(c(2,1,1),byrow=TRUE,ncol=3)

n.iter <- 1000
N <- 100
history.option <- rep(NA,n.iter)
history.value <- rep(NA,n.iter)

for(i in 1:n.iter){
	yamikumo <- FALSE
	if(i ==1 ){
		yamikumo <- TRUE
	}else{
		if(length(which(history.option[1:(i-1)]==1)) < 2 || length(which(history.option[1:(i-1)]==2)) <2){
			yamikumo <- TRUE
		}
	}
	if(yamikumo){
		s <- sample(1:2,1)
		history.option[i] <- s
		tmp.row <- sample(1:length(true.prob[[s]][,1]),1,prob=true.prob[[s]][,3])
		history.value[i] <- rnorm(1,true.prob[[s]][tmp.row,1],true.prob[[s]][tmp.row,2])
	}else{
		ones <- which(history.option ==1)
		twos <- which(history.option ==2)
		fr <- min(history.value[1:(i-1)])
		to <- max(history.value[1:(i-1)])
		fr2 <- fr-(to-fr)/2
		to2 <- to+(to-fr)/2
		d1 <- density(history.value[ones],from = fr2,to = to2)
		d2 <- density(history.value[twos],from = fr2,to = to2)
		ylim <- c(0,max(d1$y,d2$y))
		par(mfcol=c(1,1))
		plot(d1,ylim = ylim)
		par(new=TRUE)
		plot(d2,col=2,ylim=ylim)
		R1 <- sample(ones,N,replace=TRUE)
		R2 <- sample(twos,N,replace=TRUE)
		r1 <- rep(NA,N)
		r2 <- rep(NA,N)
		for(j in 1:N){
			r1[j] <- rnorm(1,history.value[R1[j]],d1$bw)
			r2[j] <- rnorm(1,history.value[R2[j]],d2$bw)
		}
		p <- length(which(r1>r2))/N
		if(runif(1)<p){
			s <- 1
		}else{
			s <- 2
		}
		history.option[i] <- s
		tmp.row <- sample(1:length(true.prob[[s]][,1]),1,prob=true.prob[[s]][,3])
		history.value[i] <- rnorm(1,true.prob[[s]][tmp.row,1],true.prob[[s]][tmp.row,2])

	}
	
}

plot(history.value, col = history.option)

多次元ベクトルデータの超球面表現と密度推定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>

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

data<-read.table("data.txt",sep="\t") #多次元ベクトルデータの読み込み(たとえば400人分、数百マーカー分)
Cdata<-cov(data,use="pairwise") # covariance matrix
Cdatacor<-cov2cor(Cdata) # correlation matrixへの変換(ただしエラーが・・・)??
answer<-princoRY(Cdata,thres=0,noPlotAxis=3,outfile="testout.txt") #princoRYの中でも『標準化』しているので、Cdataのままで与える
data4persp<-kde2d(answer$Axis.1, answer$Axis.2) #第1軸、第2軸でpersp 関数用のデータを作る。また、data4perspは密度情報
persp(data4persp) #添付図

多次元ベクトルデータの超球面表現とその先

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))
}
      • 参考の参考、Rのdensity関数
> getS3method("density","default")
function (x, bw = "nrd0", adjust = 1, kernel = c("gaussian", 
    "epanechnikov", "rectangular", "triangular", "biweight", 
    "cosine", "optcosine"), weights = NULL, window = kernel, 
    width, give.Rkern = FALSE, n = 512, from, to, cut = 3, na.rm = FALSE, 
    ...) 
{
    if (length(list(...)) > 0) 
        warning("non-matched further arguments are disregarded")
    if (!missing(window) && missing(kernel)) 
        kernel <- window
    kernel <- match.arg(kernel)
    if (give.Rkern) 
        return(switch(kernel, gaussian = 1/(2 * sqrt(pi)), rectangular = sqrt(3)/6, 
            triangular = sqrt(6)/9, epanechnikov = 3/(5 * sqrt(5)), 
            biweight = 5 * sqrt(7)/49, cosine = 3/4 * sqrt(1/3 - 
                2/pi^2), optcosine = sqrt(1 - 8/pi^2) * pi^2/16))
    if (!is.numeric(x)) 
        stop("argument 'x' must be numeric")
    name <- deparse(substitute(x))
    x <- as.vector(x)
    x.na <- is.na(x)
    if (any(x.na)) {
        if (na.rm) 
            x <- x[!x.na]
        else stop("'x' contains missing values")
    }
    N <- nx <- length(x)
    x.finite <- is.finite(x)
    if (any(!x.finite)) {
        x <- x[x.finite]
        nx <- length(x)
    }
    if (is.null(weights)) {
        weights <- rep.int(1/nx, nx)
        totMass <- nx/N
    }
    else {
        if (length(weights) != N) 
            stop("'x' and 'weights' have unequal length")
        if (!all(is.finite(weights))) 
            stop("'weights' must all be finite")
        if (any(weights < 0)) 
            stop("'weights' must not be negative")
        wsum <- sum(weights)
        if (any(!x.finite)) {
            weights <- weights[x.finite]
            totMass <- sum(weights)/wsum
        }
        else totMass <- 1
        if (!isTRUE(all.equal(1, wsum))) 
            warning("sum(weights) != 1  -- will not get true density")
    }
    n.user <- n
    n <- max(n, 512)
    if (n > 512) 
        n <- 2^ceiling(log2(n))
    if (missing(bw) && !missing(width)) {
        if (is.numeric(width)) {
            fac <- switch(kernel, gaussian = 4, rectangular = 2 * 
                sqrt(3), triangular = 2 * sqrt(6), epanechnikov = 2 * 
                sqrt(5), biweight = 2 * sqrt(7), cosine = 2/sqrt(1/3 - 
                2/pi^2), optcosine = 2/sqrt(1 - 8/pi^2))
            bw <- width/fac
        }
        if (is.character(width)) 
            bw <- width
    }
    if (is.character(bw)) {
        if (nx < 2) 
            stop("need at least 2 points to select a bandwidth automatically")
        bw <- switch(tolower(bw), nrd0 = bw.nrd0(x), nrd = bw.nrd(x), 
            ucv = bw.ucv(x), bcv = bw.bcv(x), sj = , "sj-ste" = bw.SJ(x, 
                method = "ste"), "sj-dpi" = bw.SJ(x, method = "dpi"), 
            stop("unknown bandwidth rule"))
    }
    if (!is.finite(bw)) 
        stop("non-finite 'bw'")
    bw <- adjust * bw
    if (bw <= 0) 
        stop("'bw' is not positive.")
    if (missing(from)) 
        from <- min(x) - cut * bw
    if (missing(to)) 
        to <- max(x) + cut * bw
    if (!is.finite(from)) 
        stop("non-finite 'from'")
    if (!is.finite(to)) 
        stop("non-finite 'to'")
    lo <- from - 4 * bw
    up <- to + 4 * bw
    y <- .C("massdist", x = as.double(x), xmass = as.double(weights), 
        nx = nx, xlo = as.double(lo), xhi = as.double(up), y = double(2 * 
            n), ny = as.integer(n), PACKAGE = "stats")$y * totMass
    kords <- seq.int(0, 2 * (up - lo), length.out = 2 * n)
    kords[(n + 2):(2 * n)] <- -kords[n:2]
    kords <- switch(kernel, gaussian = dnorm(kords, sd = bw), 
        rectangular = {
            a <- bw * sqrt(3)
            ifelse(abs(kords) < a, 0.5/a, 0)
        }, triangular = {
            a <- bw * sqrt(6)
            ax <- abs(kords)
            ifelse(ax < a, (1 - ax/a)/a, 0)
        }, epanechnikov = {
            a <- bw * sqrt(5)
            ax <- abs(kords)
            ifelse(ax < a, 3/4 * (1 - (ax/a)^2)/a, 0)
        }, biweight = {
            a <- bw * sqrt(7)
            ax <- abs(kords)
            ifelse(ax < a, 15/16 * (1 - (ax/a)^2)^2/a, 0)
        }, cosine = {
            a <- bw/sqrt(1/3 - 2/pi^2)
            ifelse(abs(kords) < a, (1 + cos(pi * kords/a))/(2 * 
                a), 0)
        }, optcosine = {
            a <- bw/sqrt(1 - 8/pi^2)
            ifelse(abs(kords) < a, pi/4 * cos(pi * kords/(2 * 
                a))/a, 0)
        })
    kords <- fft(fft(y) * Conj(fft(kords)), inverse = TRUE)
    kords <- pmax.int(0, Re(kords)[1:n]/length(y))
    xords <- seq.int(lo, up, length.out = n)
    x <- seq.int(from, to, length.out = n.user)
    structure(list(x = x, y = approx(xords, kords, x)$y, bw = bw, 
        n = N, call = match.call(), data.name = name, has.na = FALSE), 
        class = "density")
}
<environment: namespace:stats>
>