быстрый способ вычисления ортогонального расстояния от точки до y = x в R

У меня есть несколько точек, которые лежат вокруг y=x (см. Примеры ниже), и я надеюсь вычислить ортогональное расстояние каждой точки до этой y=x. Предположим, что точка имеет координаты (a,b), тогда легко увидеть, что спроецированная точка на y=x имеет координаты ((a+b)/2, (a+b)/2). Я использую следующие собственные коды для вычислений, но думаю, мне нужен более быстрый без for циклов. Большое спасибо!

set.seed(999)
n=50
typ.ord = seq(-2,3, length=n)   # x-axis
#
good.ord = sort(c(rnorm(n/2, typ.ord[1:n/2]+1,0.1),rnorm(n/2,typ.ord[(n/2+1):n]-0.5,0.1)))
y.min = min(good.ord)
y.max = max(good.ord)
#
plot(typ.ord, good.ord, col="green", ylim=c(y.min, y.max))
abline(0,1, col="blue")
# 
# a = typ.ord
# b = good.ord
cal.orth.dist = function(n, typ.ord, good.ord){
  good.mid.pts = (typ.ord + good.ord)/2
  orth.dist = numeric(n)
  for (i in 1:n){
    num.mat = rbind(rep(good.mid.pts[i],2), c(typ.ord[i], good.ord[i]))
    orth.dist[i] = dist(num.mat)
  }
  return(orth.dist)
}
good.dist = cal.orth.dist(50, typ.ord, good.ord)
sum(good.dist)

r distance orthogonal
person alittleboy    schedule 25.02.2013    source источник

Ответы (2)


arrow_upward
4
arrow_downward

Так же легко, как

good.dist <- sqrt((good.ord - typ.ord)^2 / 2)

Все сводится к вычислению расстояния между точкой и линией. В случае 2D y = x это становится особенно легко (попробуйте сами).

person QkuCeHBH    schedule 25.02.2013
comment
(+1) это должно быть достаточно быстро! :) - person Arun; 25.02.2013

arrow_upward
4
arrow_downward

В более общем случае (распространяясь на другие строки, возможно, более чем в двумерном пространстве), вы можете использовать следующее. Он работает путем построения матрицы проекции P из подпространства (здесь вектор A), на который вы хотите спроецировать точки x. Вычитание спроецированного компонента из точек оставляет ортогональный компонент, для которого легко вычислить расстояния.

x <- cbind(typ.ord, good.ord)          # Points to be projected
A <- c(1,1)                            # Subspace to project onto
P <- A %*% solve(t(A) %*% A) %*% t(A)  # Projection matrix P_A = A (A^T A)^-1 A^T
dists <- sqrt(rowSums(x - x %*% P)^2)  # Lengths of orthogonal residuals
person Josh O'Brien    schedule 25.02.2013
comment
Действительно мило! Но, может быть, было бы полезно проиллюстрировать, как это работает с другой линией в 2D-пространстве (другой наклон и ненулевое начало)? Не слишком уверен, как работать (то есть определять) это подпространство для проецирования на ... - person ztl; 12.01.2017
comment
@ztl Рад, что вы нашли это полезным! Для линий, проходящих через начало координат, подпространство будет определяться координатами любой точки на линии (кроме самой начала координат). Так, например, если уравнение вашей линии y = 2x, вы можете заменить c(1,1) в приведенном выше примере на c(.5, 1). Для линии, не проходящей через начало координат (то есть аффинного подпространства), вам необходимо применить алгоритм, описанный здесь. - person Josh O'Brien; 12.01.2017