Rコード最適化例:空間点パターンのメディアン Rコード最適化のコツと実例集 に戻る
Q&A コーナ投稿記事より
空間メジアン(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