2
$\begingroup$

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?

$\endgroup$
2
  • 1
    $\begingroup$ I have deleted my answer, as it turns out it was wrong! I'll redo the calculations and let you know. Sorry. $\endgroup$ – Creosote Sep 20 '15 at 11:42
  • $\begingroup$ Now redone. Enjoy. $\endgroup$ – Creosote Sep 20 '15 at 13:54
1
$\begingroup$

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.

$\endgroup$
8
  • $\begingroup$ I appreciate your answer but I am not very familiar with R but I know python. It would be absolutely great if you could post your example by python too, or write just the math if you think it is sufficient. $\endgroup$ – Dalek Sep 19 '15 at 21:31
  • $\begingroup$ are the red lines representing the distances has been computed as perpendicular distances? $\endgroup$ – Dalek Sep 20 '15 at 9:47
  • $\begingroup$ They're projections onto the "best" point on the black line. If the sd_x and sd_y values are equal for a data point, that projection should be perpendicular to the black line, or I messed up my code. $\endgroup$ – Creosote Sep 20 '15 at 9:53
  • $\begingroup$ it doesn't seem the black points are the perpendicular projections of the red points on the line. $\endgroup$ – Dalek Sep 20 '15 at 10:13
  • $\begingroup$ I could've put more comments in the code ... The black points are the true (but unseen) $x$ and $y$ values, deliberately chosen to be perfectly linear. The red circles are the observed versions of $x$ and $y$, i.e. with added "uncertainty" (Normal errors). The red lines drop down to the most-probable point on the black line. If you changed the code to sd_x=rep(2,n); sd_y=rep(0.2,n) in the relevant place, you'll see that the projections are obviously better than perpendicular would have been. $\endgroup$ – Creosote Sep 20 '15 at 10:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.