## Geometry picture
res <- 99
x <- y <- seq(0, 1, len=res)
Z0 <- Z1 <- Z2 <- matrix(NA, nrow=res, ncol=res)
mu0 <- c(0.3, 0.6)
Omega0 <- matrix(c(1,0.5, 0.5, 1), nrow=2)
mu1 <- c(0.7, 0.6)
Omega1 <- matrix(c(1, -0.5, -0.5, 1), nrow=2)
for (i in 1:res) {
  for (j in 1:res) {
    xx <- c(x[i], y[j])
    Z0[i,j] <- -crossprod(xx-mu0, Omega0)%*%(xx-mu0)
    Z1[i,j] <- -crossprod(xx-mu1, Omega1)%*%(xx-mu1)
    Omega2 <- Omega0+Omega1
    mu2 <- solve(Omega2)%*%(Omega0%*%mu0 + Omega1%*%mu1)
    Z2[i,j] <- -crossprod(xx-mu2, Omega2)%*%(xx-mu2)
  }
}
w <- seq(0, 1, len=99)
M <- matrix(NA, nrow=99, ncol=2)
for (i in 1:99) {
  Omega2 <- w[i]*Omega0+(1-w[i])*Omega1
  M[i,] <- solve(Omega2)%*%(w[i]*Omega0%*%mu0 + (1-w[i])*Omega1%*%mu1)
}
contour(x, y, Z0, drawlabels=FALSE, levels=c(-1e-3, -1e-2, -5e-2), col=pal(3)[1], xlab=expression(beta[1]), ylab=expression(beta[2]))
contour(x, y, Z1, drawlabels=FALSE, levels=c(-1e-3, -1e-2, -5e-2), col=pal(3)[3], add=TRUE)
contour(x, y, Z2, drawlabels=FALSE, levels=c(-1e-3, -1e-2, -5e-2), col=pal(3)[2], add=TRUE)
lines(M[,1], M[,2], lwd=3)