title: Reproduce Fisher linear discriminant figure
body: <p>Many books illustrate the idea of Fisher linear discriminant analysis using the following figure (this particular is from <a href="http://research.microsoft.com/en-us/um/people/cmbishop/prml/" rel="nofollow noreferrer">Pattern Recognition and Machine Learning</a>, p. 188)</p>
  
  <p><a href="https://i.stack.imgur.com/4afF3.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/4afF3.jpg" alt="enter image description here"></a></p>
  
  <p>I wonder how to reproduce this figure in R (or in any other language). Pasted below is my initial effort in R. I simulate two groups of data and draw linear discriminant using <code>abline()</code> function. Any suggestions are welcome.</p>
  
  <pre><code>set.seed(2014)
  library(MASS)
  library(DiscriMiner) # For scatter matrices
  
  # Simulate bivariate normal distribution with 2 classes
  mu1 &lt;- c(2, -4)
  mu2 &lt;- c(2, 6)
  rho &lt;- 0.8
  s1 &lt;- 1
  s2 &lt;- 3
  Sigma &lt;- matrix(c(s1^2, rho * s1 * s2, rho * s1 * s2, s2^2), byrow = TRUE, nrow = 2)
  n &lt;- 50
  X1 &lt;- mvrnorm(n, mu = mu1, Sigma = Sigma)
  X2 &lt;- mvrnorm(n, mu = mu2, Sigma = Sigma)
  y &lt;- rep(c(0, 1), each = n)
  X &lt;- rbind(x1 = X1, x2 = X2)
  X &lt;- scale(X)
  
  # Scatter matrices
  B &lt;- betweenCov(variables = X, group = y)
  W &lt;- withinCov(variables = X, group = y)
  
  # Eigenvectors
  ev &lt;- eigen(solve(W) %*% B)$vectors
  slope &lt;- - ev[1,1] / ev[2,1]
  intercept &lt;- ev[2,1]
  
  par(pty = "s")
  plot(X, col = y + 1, pch = 16)
  abline(a = slope, b = intercept, lwd = 2, lty = 2)
  </code></pre>
  
  <p><strong>MY (UNFINISHED) WORK</strong></p>
  
  <p>I pasted my current solution below. The main question is how to rotate (and move) the density plot according to decision boundary. Any suggestions are still welcome.</p>
  
  <p><a href="https://i.stack.imgur.com/4Z4DB.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/4Z4DB.png" alt="enter image description here"></a></p>
  
  <pre><code>require(ggplot2)
  library(grid)
  library(MASS)
  
  # Simulation parameters
  mu1 &lt;- c(5, -9)
  mu2 &lt;- c(4, 9)
  rho &lt;- 0.5
  s1 &lt;- 1
  s2 &lt;- 3
  Sigma &lt;- matrix(c(s1^2, rho * s1 * s2, rho * s1 * s2, s2^2), byrow = TRUE, nrow = 2)
  n &lt;- 50
  # Multivariate normal sampling
  X1 &lt;- mvrnorm(n, mu = mu1, Sigma = Sigma)
  X2 &lt;- mvrnorm(n, mu = mu2, Sigma = Sigma)
  # Combine into data frame
  y &lt;- rep(c(0, 1), each = n)
  X &lt;- rbind(x1 = X1, x2 = X2)
  X &lt;- scale(X)
  X &lt;- data.frame(X, class = y)
  
  # Apply lda()
  m1 &lt;- lda(class ~ X1 + X2, data = X)
  m1.pred &lt;- predict(m1)
  # Compute intercept and slope for abline
  gmean &lt;- m1$prior %*% m1$means
  const &lt;- as.numeric(gmean %*% m1$scaling)
  z &lt;- as.matrix(X[, 1:2]) %*% m1$scaling - const
  slope &lt;- - m1$scaling[1] / m1$scaling[2]
  intercept &lt;- const / m1$scaling[2]
  
  # Projected values
  LD &lt;- data.frame(predict(m1)$x, class = y)
  
  # Scatterplot
  p1 &lt;- ggplot(X, aes(X1, X2, color=as.factor(class))) + 
    geom_point() +
    theme_bw() +
    theme(legend.position = "none") +
    scale_x_continuous(limits=c(-5, 5)) + 
    scale_y_continuous(limits=c(-5, 5)) +
    geom_abline(intecept = intercept, slope = slope)
  
  # Density plot 
  p2 &lt;- ggplot(LD, aes(x = LD1)) +
    geom_density(aes(fill = as.factor(class), y = ..scaled..)) +
    theme_bw() +
    theme(legend.position = "none")
  
  grid.newpage()
  print(p1)
  vp &lt;- viewport(width = .7, height = 0.6, x = 0.5, y = 0.3, just = c("centre"))
  pushViewport(vp)
  print(p2, vp = vp)
  </code></pre>