Including errors in perpendicular off-set determination of points from a line

Question

I have a set of lines and some data points which have errors for the x and y components I would like to see which one of these lines have the minimum perpendicular offset from the points. so I can find the perpendicular distances without including the errors as following $$R=\sum_{i=1}^{N}\frac{|y_i-(a+bx_i)|}{\sqrt{1+b^2}}$$ and find the minimum value of $R$ in order to get the best line. I am wondering how I can include $x_i$ and $y_i$ uncertainties?

Answer 1

This is a completely reworked answer. Apologies for previous errors.

You have a data point $(x,y)$, which is a "fuzzy" look at some unknown $(x^\ast,y^\ast)$ such that $x = x^\ast+\textrm{Normal}(0,\sigma_x^2)$ and $y=y^\ast+\textrm{Normal}(0,\sigma_y^2)$. You want to see if $(x,y)$ look like they came from the line $Y=a+bX$, taking into account the uncertainty $(\sigma_x,\sigma_y)$.

With that set-up, the best-matching $X$ on the line $Y=a+bX$ is $$ \dfrac{b(y-a)\sigma_x^2+x\sigma_y^2}{b^2\sigma_x^2+\sigma_y^2} $$ and the log of the probability density at that point is $$ \textrm{constant} - \dfrac{1}{2}\frac{(y-a-bx)^2}{b^2\sigma_x^2+\sigma_y^2} $$

The R code shown below finds $(a,b)$ maximising that probability density. It's drawn as a red line. For comparison, the ordinary least-squares fit is shown as a blue line. True $(x^\ast,y^\ast)$ are in black -- these are what you'd like to find with your line -- and observed $(x,y)$ are in red. The red and blue lines are often about equally close to the true data. Whenever they are vastly different, though, it's always blue (ordinary least squares) that's worse, which I take as confirmation that the "red method" is doing the right thing with those $(\sigma_x,\sigma_y)$ uncertainties.

# Score a line Y = a + b*X against a set of data points {(x,y)}
# each subject to N(0,sd_x) and N(0,sd_y) noise.
# (lower scores are better)
#
line_score = function(x,y,sd_x,sd_y,a,b) {
  sum( (y-a-b*x)**2 / (b**2*sd_x**2+sd_y**2) )
}

# Find the X on Y = a + b*X that maximises the probability of a
# particular point (x,y) with N(0,sd_x) and N(0,sd_y) uncertainty.
#
best_X = function(x,y,sd_x,sd_y,a,b) {
  (b*(y-a)*sd_x**2+x*sd_y**2) / (b**2*sd_x**2+sd_y**2)
}

# Fit a line through data points {(x,y)} each subject to N(0,sd_x)
# and N(0,sd_y) noise.
#
b_error = function(x,y,sd_x,sd_y,bb) {
  zz = (bb**2*sd_x**2+sd_y**2)
  aa = sum((y-bb*x)/zz) / sum(1/zz)
  sum( (y-aa-bb*x)**2 / zz )
}
fit = function(x,y,sd_x,sd_y) {
  b1 = 10
  while (b_error(x,y,sd_x,sd_y,b1+1) < b_error(x,y,sd_x,sd_y,b1)) b1=b1+1
  b1 = b1+1
  b0 = b1
  while (b_error(x,y,sd_x,sd_y,b0-1) < b_error(x,y,sd_x,sd_y,b0)) b0=b0-1
  b0 = b0-1
  while (abs(b0/b1-1) > 1e-3) {
    em = b_error(x,y,sd_x,sd_y,(b0+b1)/2)
    bu = (3*b0+b1)/4; eu = b_error(x,y,sd_x,sd_y,bu)
    bv = (b0+3*b1)/4; ev = b_error(x,y,sd_x,sd_y,bv)
    if (eu>=em) { b0=bu }
    if (ev>em) { b1=bv }
  }
  bb = (b0+b1)/2
  zz = (bb**2*sd_x**2+sd_y**2)
  aa = sum((y-bb*x)/zz) / sum(1/zz)
  list(a=aa,b=bb)
}
fit(x,y,sd_x,sd_y)

# Create some sample data.
#
make_data = function(a,b,n) {
  tx=seq(1:n)
  ty=a+b*tx
  x=tx+rnorm(n,0,sd_x)
  y=ty+rnorm(n,0,sd_y)
  list(x=x,y=y,tx=tx,ty=ty)
}

# Plot the data.
#
plot_data = function(x,y,tx,ty,a,b) {
  tmp = lm(y~x)
  lma = tmp$coefficients[1]
      lmb = tmp$coefficients[2]
  rrls = line_score(x,y,sd_x,sd_y,lma,lmb)
  tmp = fit(x,y,sd_x,sd_y)
  fla = tmp$a
      flb = tmp$b
  fls = line_score(x,y,sd_x,sd_y,fla,flb)
  lim=range(c(x,y,tx,ty))
  plot(
    tx,ty,pch=20,xlim=lim,ylim=lim,xlab="X",ylab="Y",
    main=paste("line_score(R**2 line) =",rrls,"\nline_score(fitted line) =",fls)
  )
  points(x,y,col="red")
  show_lines = function(a,b,col) {
    abline(a,b,col=col)
    for (i in 1:n) {
      lx = best_X(x[i],y[i],sd_x[i],sd_y[i],a,b)
      ly = a + b*lx
      lines(c(x[i],lx),c(y[i],ly),lty=3,col=col)
    }
  }
  show_lines(fla,flb,"red")
  show_lines(lma,lmb,"blue")
}

# Parameters for the experiment.
#
a=2; b=1; n=10
sd_x=runif(n)*3; sd_y=runif(n)/2

# Do it.
#
data = make_data(a,b,n)
plot_data(data$x,data$y,data$tx,data$ty,a,b)

---

Supplementary material ...

With $\textrm{Normal}(0,\sigma_x^2)$ uncertainty around $x$, a point $X$ near $x$ has log probability density equal to a constant (independent of $X$) minus $(X-x)^2/(2\sigma_x^2)$. Similarly for $y$ and $Y$. So the best $X$ on the line $Y=a+bX$ minimises $(\frac{X-x}{\sigma_x})^2+(\frac{a+bX-y}{\sigma_y})^2$. Differentiating with respect to $X$ and equating to zero gives the best-matching $X$ shown in the main part of this answer. Note, by the way, that if $y=a+bx$ (i.e. the point is on the line) then the best-matching $X$ is $x$ as it should be.