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

SIZE(20){COLOR(magenta){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

タイムスタンプを変更しない

SIZE(20){COLOR(magenta){Rコード最適化例：空間点パターンのメディアン}}  [[Rコード最適化のコツと実例集]] に戻る~

テキスト整形のルールを表示する

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

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