ryamadaのコンピュータ・数学メモ

2014-01-26

傾きを単調にする

$X=(x_1,...,x_n),Y=(y_1,...,y_n)$ という点列がある。X,Yともに単調増であるとする。
ここに傾きが単調減な線を引こうと思ったらIsotonic regressionすることで可能
Rでは

library(quadprog)
# データを適当に作る
n <- 100
x <- sort(runif(n))
y <- cumsum(sort(runif(n),decreasing=TRUE))
plot(x,y)

傾きが単調減であるとは $\frac{y_{i+2}-y_{i+1}}{x_{i+2}-x_{i+1}} \le \frac{y_{i+1}-y_i}{x_{i+1}-x_i}$
この式を変形して[tex:(x_{i+1}-x_i)y_{i+2}-*1y_{i+1}-(x_{i+2}-x_{i+1})y_i \le 0]

# 傾きが単調な線を引く
# 各点の重みを平等にするための行列
Dmat <- diag(rep(1,n))
# 制約行列は制約数n-2行、レコード数n列
diff.x <- diff(x)
Amat <- matrix(0,n-2,n)
for(i in 1:(n-2)){
	Amat[i,i] <- -diff.x[i+1]
	Amat[i,i+1] <- (diff.x[i]+diff.x[i+1])
	Amat[i,i+2] <- -diff.x[i]
}
# Amatは制約不等式の係数。その不等式の右辺の値を納めたベクトル
bvec <- rep(0,n-2)
# meqは制約のうちの等式の数(今回はない)
# 解く
out <- solve.QP(Dmat,y,t(Amat),bvec=bvec,meq=0)
# プロットする
plot(x,y,pch=20,cex=0.8)
points(x,out[[1]],type="l",col=2)

*1:x_{i+2}-x_{i+1})+(x_{i+1}-x_i