ryamadaの遺伝学・遺伝統計学メモ このページをアンテナに追加 RSSフィード

数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

2014-03-27 カーネル関数とオーダー

[][]カーネル推定のオーダー

  • こちらで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であるような関数であると言うことにする
    • 対称関数ではqには偶数のみが適合し、2-th orderカーネル,4-th orderカーネル…となる
  • ちなみに、非負の値しかとらないカーネル関数では、¥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}}

f:id:ryamada22:20140325090247j:image

  • じゃあ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)

f:id:ryamada22:20140325091415j:image

  • さらに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)

f:id:ryamada22:20140325093446j:image

  • 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}

f:id:ryamada22:20140325093447j:image

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")

2014-03-26 カーネル推定と特性関数とモーメントとフーリエ変換・級数

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

  • こちらで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での式・値になる
  • この積率母関数は、確率密度関数そのものと同様に、確率変数をきちんと定義づける関数であるが、そのような性質を持つ関数に特性関数がある
    • 特性関数は¥phi_X(t) = E¥[e^{itX}¥]のように指数が複素数になったものであって、積率母関数と同様にモーメントをもたらす
      • この表現自体が「フーリエ変換」と見える(参考)
    • E(X^n) = i^{-n} ¥phi_X^{(n)}(0) = i^{-n} ¥[¥frac{d^n ¥phi_X(t)}{dt^n}¥]である
  • 確率密度関数にある特性関数はカーネル関数にも定義できる(多分)。そうすると、モーメントの期待値が特性関数の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の構図に持ち込まれる

2013-03-01 カーネル密度推定の上下限

ryamada222013-03-01

[][][]折り返し

f:id:ryamada22:20130301045334j:image

x <- seq(from=-10,to=10,by=0.01)
m <-1.3
y <- dnorm(x,m)
posi <- which(x>=0)
neg <- which(x<=0)
yposi <- y[posi]
yneg <- y[neg[length(neg):1]]
y1 <- yposi+yneg

plot(y1,type="l")

2013-02-28 量的観察に基づく決断

ryamada222013-02-28

[][][]量的データで決断

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

f:id:ryamada22:20130227161935j:image

  • 書き散らしソース
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)


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

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