gg.mixEM <- function(EM) {
  require(ggplot2)
  x       <- with(EM,seq(min(x),max(x),len=1000))
  pars    <- with(EM,data.frame(comp=colnames(posterior), mu, sigma,lambda))
  em.df   <- data.frame(x=rep(x,each=nrow(pars)),pars)
  em.df$y <- with(em.df,lambda*dnorm(x,mean=mu,sd=sigma))
  ggplot(data.frame(x=EM$x),aes(x,y=..density..)) +
      geom_histogram(fill=NA,color="grey")+
      geom_path(data=em.df,aes(x,y,color=comp),alpha=1)+
      scale_color_discrete("Component\nMeans",labels=format(em.df$mu,digits=3))+
      theme_bw()
}

intersect <- function(m1, sd1, m2, sd2, p1, p2){
  B <- (m1/sd1^2 - m2/sd2^2)
  A <- 0.5*(1/sd2^2 - 1/sd1^2)
  C <- 0.5*(m2^2/sd2^2 - m1^2/sd1^2) - log((sd1/sd2)*(p2/p1))

  if (A!=0){
      (-B + c(1,-1)*sqrt(B^2 - 4*A*C))/(2*A)
  } else {-C/B}
}

library(mixtools)

tauvals=pbpastex()$V1

tauvals.nona = as.numeric(tauvals[!is.na(as.numeric(tauvals))])
mixmodel = normalmixEM(tauvals.nona)

intersect_point = intersect(mixmodel$mu[1],mixmodel$sigma[1],mixmodel$mu[2],mixmodel$sigma[2],mixmodel$lambda[1],mixmodel$lambda[2])[2]

gg.mixEM(mixmodel)
  +geom_vline(xintercept=intersect_point)
  +ggtitle(paste("Cutoff",intersect_point))
  +geom_path(data=data.frame(x=tauvals.nona),aes(x=x),stat="density",linetype="dotted")