mixband.R

mixband <- function(y){
    ## y1 <- scale(y)
    y1 <- y
    m1 <- tryCatch(normalmixEM(y1),error=function(e) return(NULL))
    if(is.null(m1)) return(c(l=0,pr=0, p2=0,isBandy2=0))
    lambda <- m1$lambda
    means <- m1$mu
    clusters <- apply(m1$posterior,1,function(r) if(r[1]>r[2]) 1 else 2)
    if(any( table(clusters)<=2) || length(table(clusters))==1 ) return(c(l=max(lambda)/min(lambda),pr=0,p2=0,isBandy2=0))
    prsep <- min(as.numeric(tapply(apply(m1$posterior,1,function(r) max(r)-min(r)),clusters, mean)))
    icdf <- qnorm(seq_along(y)/(length(y)+1))
    (t2 <- t.test( icdf[ clusters==1], icdf[clusters==2] ,alternative="two.sided",var.equal=FALSE))
    isBandy2 <- 1*(1*(t2$p.value>0.01) && max(lambda)/min(lambda)<=5 && prsep>=0.90)
    c(l=max(lambda)/min(lambda), pr=prsep,p2=1*(t2$p.value>0.01), isBandy2=isBandy2)
}