Rコード最適化例：空間点パターンのメディアン

Rコード最適化例：空間点パターンのメディアン Rコード最適化のコツと実例集 に戻る

Ｑ＆Ａ コーナ投稿記事より
空間メジアン(spatial median)と2変量オヤ・メジアン(bivariate Oja median)を計算する。 T(θ) = Σ|xi-θ|とするとき T(θ) を最小にするθが空間メジアン、 T(θ)=ΣA(xi,xj,θ)とするとき T(θ) を最小にするθが2変量オヤ・メジアン (A(a,b,c)は3点a,b,cを頂点とする三角形の面積)。

## 空間メディアンを optim() で求める
> X <- runif(1000)  # テスト用の1000個の点の x 座標
> Y <- runif(1000)  # テスト用の1000個の点の y 座標
## 任意の点 T=(T[1], T[2]) から点 (X[i],Y[i]) までの距離和関数
dsum <- function (T) {  # 素朴な定義
          s <- 0
          for (i in 1:length(X)) s <- s + sqrt((T[1]-X[i])^2+(T[2]-Y[i])^2)
          s }
> system.time(res <- optim(c(mean(X),mean(Y)), dsum))   # 重心を初期値にとり Nelder-Mead 法で最小化
[1] 0.62 0.00 1.20 0.00 0.00                            # 関数評価回数41回
> dsum2 <- function(T) { sum(sqrt((T[1]-X)^2+(T[2]-Y)^2)) } # ベクトル化演算で高速化
> system.time(res2 <- optim(c(mean(X),mean(Y)), dsum2)) # 重心を初期値にとり Nelder-Mead 法で最小化
[1] 0.03 0.00 0.07 0.00 0.00

## 2変量オヤ・メジアンを optim() で求める 
## 三点 (X[i],Y[i]),(X[j],Y[j]), (T[1], T[2]) の作る三角形の面積は 
##   abs((X[i]-T[i])*(Y[j]-T[j]) - (Y[i]-T[i])*(X[j]-T[j]))/2 に注意する
## 三角形面積和を計算する素朴な関数
asum <- function (T) {
          asum <- 0
          for (i in 1:(length(X)-1)) # 全ての 1 <= i < j <= length(X) に対する和
            for (j in (i+1):length(X)) 
              asum <- asum + abs((X[i]-T[1])*(Y[j]-T[2]) - (Y[i]-T[1])*(X[j]-T[2]))
          asum/2 }   # 毎回割らずに、最後に割っていることを注意
## ベクトルの外積演算を使い高速化
## (X-T[1])%o%(Y-T[2]) は (i,j) 成分が (X[i]-T[1])*(Y[j]-T[2]) である行列
> asum2 <- function (T) {Z <- (X-T[1])%o%(Y-T[2]); sum(abs((Z-t(Z))))/4}
> asum2(c(0,0))  # 検査
[1] 665.7858
> X <- runif(100) # テスト用の100個の点の x 座標
> Y <- runif(100) # テスト用の100個の点の y 座標
> system.time(res <- optim(c(mean(X),mean(Y)), asum)) # 重心を初期値にとり Nelder-Mead 法で最小化
[1]  5.6  0.0 11.1  0.0  0.0  # 関数評価回数 47 回
> system.time(res2 <- optim(c(mean(X),mean(Y)), asum2)) # 重心を初期値にとり Nelder-Mead 法で最小化
[1] 0.09 0.00 0.14 0.00 0.00

>  X <- runif(200) # 点の数を200にする
>  Y <- runif(200)
> system.time(res2 <- optim(c(mean(X),mean(Y)), asum2)) # 重心を初期値にとり Nelder-Mead 法で最小化
[1] 0.43 0.03 0.86 0.00 0.00  # 関数評価回数47回
> system.time(res <- optim(c(mean(X),mean(Y)), asum)) # 重心を初期値にとり Nelder-Mead 法で最小化
[1] 21.11  0.29 43.64  0.00  0.00

最新の20件