# R/findInt.r In MarcelaCespedes/FSPmix: Exploratory tool for fast unsupervised features selection and clustering using mixture models

```# relies on this function
# @@ Alternative method - just find the intersection between the two fitted Gaussians.
#' Find intersection(s) between two Gaussians, possibly with mixing proportions.
#' @param m1, m2: fitted means of two components
#' @param sd1, sd2: fitted standard deviations of two components
#' @param p1, p2: fitted proportions for the two components
#' @param filter {boolean}: if TRUE, try to find the "middle" & real-valued intersection.
#' There can be more than one intersection between two Gaussians,
#' We try to find the Real one that is between the two means.
#' i.e. the "middle" one.
#' I have the working .. somewhere .., I dug this out of my code.
#' Similar to this:: (but with mixing proportions)
#' https://stats.stackexchange.com/questions/311592/how-to-find-the-point-where-two-normal-distributions-intersect
#' @examples
#' library(mixtools)
#' y <- rnormmix(1e3, lambda=c(0.2, 0.8), mu=c(-0.5, 0.5), sigma=c(0.1, 0.3))
#' hist(y)
#' m <- normalmixEM(y, k=2)
#' thr <- findInt(m\$mu, m\$mu, m\$sigma, m\$sigma, m\$lambda, m\$lambda)
#' DenPlot(m, BiomThresh=thr, BiomName='y', BinWidth=0.1, BiomRange=waiver())
findInt <- function (m1, m2, sd1, sd2, p1=1, p2=1, filter=T) {
a <- 1/(2*sd1^2) - 1/(2*sd2^2)
b <- m2/(sd2^2) - m1/(sd1^2)
c <- m1^2 /(2*sd1^2) - m2^2 / (2*sd2^2) - log(sd2/sd1)
c <- m1^2 /(2*sd1^2) - m2^2 / (2*sd2^2) - log((sd2*p1)/(sd1*p2))
rr <- polyroot(c(c,b,a))

sw = 0
# find the one between the means
if (filter) {
rr <- Re(rr[sapply(Im(rr), function (x) isTRUE(all.equal(x, 0)))])
rr <- rr[rr >= min(m1, m2) & rr <= max(m1, m2)]
#rr <- rr[rr >= m1 & rr <= m2]
if (length(rr) != 1){
sw = 1
#warning(sprintf("Failed to find intersection between the means (%i candidates)",
#             length(rr))) # supress warnings, throw out switch instead
#stop(sprintf("Failed to find intersection between the means (%i candidates)", length(rr)))
}
}
return(list(rr = rr, sw = sw))
}
```
MarcelaCespedes/FSPmix documentation built on May 12, 2020, 5:49 p.m.