2013-10-07

ドロネー分割を取り入れる

昨日の記事では酔歩を使って平滑化する話
そのとき、観測点の密度が違うと平滑化後の値が密度推定にはならないことを書いた
それを修正するために、ドロネー図的に(多次元に一般化しやすくするため)各点に与えうる密度っぽい値を計算することにしたこちらにドロネー図やら、ドロネー図的体積計算やらを書いた
もちろん、それを使うと、凸凹になるので、密度っぽい値にした後で酔歩平滑化を入れることにする
まずは、１次元分布について、いろいろな例を見てみよう
単純な正規分布。図は、青が「真の分布」、黒はRのdensity()によるSJ-ste幅による推定密度、赤はドロネー体積を使わず、ただの酔歩平滑化による値、緑はドロネー体積による補正後の酔歩平滑化による値
どれも同じような結果

# 単体の体積
my.simplex.volume <- function(V,factorial=FALSE){
	n <- length(V[1,])
	K <- matrix(V[1:n,],ncol=n)
	K. <- t(K)-V[n+1,]
	if(factorial){
		return(det(K.)/factorial(n))
	}else{
		return(det(K.))
	}
	
}

# 単体の体積計算をドロネー頂点ID情報で行うための関数
my.delaunay.vol <- function(a,X){
	tmp <- matrix(X[a,],ncol=length(X[1,]))
	my.simplex.volume(tmp)
}
# Xは観測点座標行列、n,rは酔歩平滑化のためのパラメタ
my.delaunay.rw.density <- function(X,n=1000,r=0.1){
	# ドロネー化(１次元の場合は、大小隣の点と結び合う線分が指定される)
	delaunay.X <- delaunayn(X)
	# ドロネー三角形(の一般次元化)体積を計算
	# １次元の場合は線分長
	spx.vol <- apply(delaunay.X,1,my.delaunay.vol,X)
	# 各観察点の「ドロネー的体積」
	Vol <- rep(0,length(X[,1]))
	for(i in 1:length(Vol)){
		tmp <- (delaunay.X==i)
		tmp.2 <- apply(tmp,1,any)
		Vol[i] <- sum(abs(spx.vol[which(tmp.2)]))
	}
	# 「ドロネー的体積」で除することで、１観察の空間的重みに換算する
	Y.d <- 1/Vol
	# 酔歩による平滑化をこの「ドロネー的体積補正」した後の観測状態に実施する
	RW.Mx.d <- my.rw.density.mx(X,Y.d,n,r)
	# 平滑化後の各点の重み
	Y.post <- RW.Mx.d %*% Y.d
	# Volは各点のドロネー的体積、Y.dはその逆数。各点での密度に相当
	# Y.postは平滑化後の各点の重み
	# Mは平滑のための行列(指数行列処理にて酔歩平衡状態での移動行列)
	return(list(Vol=Vol,Y.d=Y.d,Y.post=Y.post,M = RW.Mx.d))
}

# １次元の分布
# 母分布を決めよう
means <- c(0,0)
sds <- c(1,1)
ratio <- c(0.4,0.6)
N.pt <- 50
X <- c(rnorm(N.pt*ratio[1],means[1],sds[1]),rnorm(N.pt*ratio[2],means[2],sds[2]))
true.Dens <- ratio[1] * dnorm(X,means[1],sds[1])+ratio[2]*dnorm(X,means[2],sds[2])
# 行列型にする
X <- matrix(X,ncol=1)
# 各点に１観測ずつ
Y <- rep(1,length(X[,1]))

# 酔歩による平滑化をこの「ドロネー的体積補正」する前とした後の観測状態に実施するn <- 1000
r <- 0.1
RW.Mx <- my.rw.density.mx(X,Y,n,r)
# 平滑化後の各点の重み
Z <- RW.Mx %*% Y

# ドロネー的体積補正法
D.out <- my.delaunay.rw.density(X,n,r)
Z.d <- D.out$Y.post

dens.X <- density(X,bw = "SJ-ste")
plot(dens.X,lwd=3)
ord <- order(X)
points(X[ord],Z[ord]/max(Z)*max(dens.X[[2]]),pch=20,col=2,type="b")
points(X[ord],Z.d[ord]/max(Z.d)*max(dens.X[[2]]),pch=20,col=3,type="b")
points(X[ord],true.Dens[ord]/max(true.Dens)*max(dens.X[[2]]),pch=20,col=4,type="l",lwd=3)

２峰性にしてみる
- density()が２峰の谷をうまく拾っているようだ

# １次元の分布
# 母分布を決めよう
means <- c(0,3)
sds <- c(1,1)
ratio <- c(0.4,0.6)
N.pt <- 200
X <- c(rnorm(N.pt*ratio[1],means[1],sds[1]),rnorm(N.pt*ratio[2],means[2],sds[2]))
true.Dens <- ratio[1] * dnorm(X,means[1],sds[1])+ratio[2]*dnorm(X,means[2],sds[2])
# 行列型にする
X <- matrix(X,ncol=1)
# 各点に１観測ずつ
Y <- rep(1,length(X[,1]))

# 酔歩による平滑化をこの「ドロネー的体積補正」する前とした後の観測状態に実施するn <- 1000
r <- 0.1
RW.Mx <- my.rw.density.mx(X,Y,n,r)
# 平滑化後の各点の重み
Z <- RW.Mx %*% Y

# ドロネー的体積補正法
D.out <- my.delaunay.rw.density(X,n,r)
Z.d <- D.out$Y.post

dens.X <- density(X,bw = "SJ-ste")
plot(dens.X,lwd=3)
ord <- order(X)
points(X[ord],Z[ord]/max(Z)*max(dens.X[[2]]),pch=20,col=2,type="b")
points(X[ord],Z.d[ord]/max(Z.d)*max(dens.X[[2]]),pch=20,col=3,type="b")
points(X[ord],true.Dens[ord]/max(true.Dens)*max(dens.X[[2]]),pch=20,col=4,type="l",lwd=3)

もっと離してみる
- 「離す」と酔歩平滑化法では、離れた２つの峰の間の酔歩拡散が起きにくくなるので、ドロネー的密度補正をしないと(赤)、山の高さが２峰で変わらない問題が起きるが、ドロネー的密度補正でそれが防げていることがわかる

# １次元の分布
# 母分布を決めよう
means <- c(0,10)
sds <- c(1,1)
ratio <- c(0.4,0.6)
N.pt <- 200
X <- c(rnorm(N.pt*ratio[1],means[1],sds[1]),rnorm(N.pt*ratio[2],means[2],sds[2]))
true.Dens <- ratio[1] * dnorm(X,means[1],sds[1])+ratio[2]*dnorm(X,means[2],sds[2])
# 行列型にする
X <- matrix(X,ncol=1)
# 各点に１観測ずつ
Y <- rep(1,length(X[,1]))

# 酔歩による平滑化をこの「ドロネー的体積補正」する前とした後の観測状態に実施するn <- 1000
r <- 0.1
RW.Mx <- my.rw.density.mx(X,Y,n,r)
# 平滑化後の各点の重み
Z <- RW.Mx %*% Y

# ドロネー的体積補正法
D.out <- my.delaunay.rw.density(X,n,r)
Z.d <- D.out$Y.post

dens.X <- density(X,bw = "SJ-ste")
plot(dens.X,lwd=3)
ord <- order(X)
points(X[ord],Z[ord]/max(Z)*max(dens.X[[2]]),pch=20,col=2,type="b")
points(X[ord],Z.d[ord]/max(Z.d)*max(dens.X[[2]]),pch=20,col=3,type="b")
points(X[ord],true.Dens[ord]/max(true.Dens)*max(dens.X[[2]]),pch=20,col=4,type="l",lwd=3)

離れた２ピークの分散に差が出てくると分布推定は難しくなってくる

# １次元の分布
# 母分布を決めよう
means <- c(0,50)
sds <- c(1,8)
ratio <- c(0.2,0.8)
N.pt <- 500
X <- c(rnorm(N.pt*ratio[1],means[1],sds[1]),rnorm(N.pt*ratio[2],means[2],sds[2]))
true.Dens <- ratio[1] * dnorm(X,means[1],sds[1])+ratio[2]*dnorm(X,means[2],sds[2])
# 行列型にする
X <- matrix(X,ncol=1)
# 各点に１観測ずつ
Y <- rep(1,length(X[,1]))

# 酔歩による平滑化をこの「ドロネー的体積補正」する前とした後の観測状態に実施するn <- 1000
r <- 0.1
RW.Mx <- my.rw.density.mx(X,Y,n,r)
# 平滑化後の各点の重み
Z <- RW.Mx %*% Y

# ドロネー的体積補正法
D.out <- my.delaunay.rw.density(X,n,r)
Z.d <- D.out$Y.post

dens.X <- density(X,bw = "SJ-ste")
#dens.X <- density(X)
plot(dens.X,lwd=3)
ord <- order(X)
points(X[ord],Z[ord]/max(Z)*max(dens.X[[2]]),pch=20,col=2,type="b")
points(X[ord],Z.d[ord]/max(Z.d)*max(dens.X[[2]]),pch=20,col=3,type="b")
points(X[ord],true.Dens[ord]/max(true.Dens)*max(dens.X[[2]]),pch=20,col=4,type="l",lwd=3)

分散の異なる分布が重なっているとき

# １次元の分布
# 母分布を決めよう
means <- c(0,2)
sds <- c(1,8)
ratio <- c(0.3,0.7)
N.pt <- 500
X <- c(rnorm(N.pt*ratio[1],means[1],sds[1]),rnorm(N.pt*ratio[2],means[2],sds[2]))
true.Dens <- ratio[1] * dnorm(X,means[1],sds[1])+ratio[2]*dnorm(X,means[2],sds[2])
# 行列型にする
X <- matrix(X,ncol=1)
# 各点に１観測ずつ
Y <- rep(1,length(X[,1]))

# 酔歩による平滑化をこの「ドロネー的体積補正」する前とした後の観測状態に実施するn <- 1000
r <- 0.1
RW.Mx <- my.rw.density.mx(X,Y,n,r)
# 平滑化後の各点の重み
Z <- RW.Mx %*% Y

# ドロネー的体積補正法
D.out <- my.delaunay.rw.density(X,n,r)
Z.d <- D.out$Y.post

dens.X <- density(X,bw = "SJ-ste")
#dens.X <- density(X)
plot(dens.X,lwd=3)
ord <- order(X)
points(X[ord],Z[ord]/max(Z)*max(dens.X[[2]]),pch=20,col=2,type="b")
points(X[ord],Z.d[ord]/max(Z.d)*max(dens.X[[2]]),pch=20,col=3,type="b")
points(X[ord],true.Dens[ord]/max(true.Dens)*max(dens.X[[2]]),pch=20,col=4,type="l",lwd=3)