## Illustration of how OLS and total least squares compare for the
## simple linear model.

## Samples size.
n = 100

for (k in 1:100) {

  ## Generate data from a simple linear model.
  X = rnorm(n)
  Y = 1 + X + rnorm(100)

  ## OLS regression.
  m = lm(Y ~ X)

  ## Total least squares regression.
  Q = cbind(X, Y)
  C = cov(Q)

  ## B is normal to the regression line.
  E = eigen(C)
  B = E$vectors[,2]

  ## D is parallel to the regression line.
  D = c(-B[2]/B[1], 1)
  
  ## Center of the data.
  M = apply(Q, 2, mean)

  ## Convert the TLS line to y = alpha + beta*x form.
  v = -M[1]/D[1]
  alpha = M[2] + v*D[2]
  beta = D[2]/D[1]
  
  ## Plot the data.
  plot(X, Y, t='p', ylab="Y", xlab="X")

  ## Plot the OLS fitted line.
  lines(X, m$fitted.values, t='l', col='red')

  ## Plot the TLS fitted line.
  lines(X, alpha + beta*X, t='l', col='blue')

  legend(x='topleft', legend=c('OLS', 'TLS'),
         col=c('red', 'blue'), lty=c(1,1))
  
  scan()
}