qnormMix rewrite fix #2

DESCRIPTION

@@ -2,7 +2,7 @@ Package: EnvStats
 Type: Package
 Title: Package for Environmental Statistics, Including US EPA Guidance
 Version: 2.4.0
-Date: 2019-01-14
+Date: 2020-10-20
 Authors@R: c(
   person("Steven P.","Millard",email="EnvStats@ProbStatInfo.com",role=c("aut")),
   person("Alexander", "Kowarik", email = "alexander.kowarik@statistik.gv.at", role = c("ctb","cre"),

R/qnormMix.R

@@ -1,45 +1,14 @@
-qnormMix <-
-function (p, mean1 = 0, sd1 = 1, mean2 = 0, sd2 = 1, p.mix = 0.5) 
-{
-    names.p <- names(p)
-    arg.mat <- cbind.no.warn(p = as.vector(p), mean1 = as.vector(mean1), 
-        sd1 = as.vector(sd1), mean2 = as.vector(mean2), sd2 = as.vector(sd2), 
-        p.mix = as.vector(p.mix))
-    na.index <- is.na.matrix(arg.mat)
-    if (all(na.index)) 
-        q <- rep(NA, nrow(arg.mat))
-    else {
-        q <- numeric(nrow(arg.mat))
-        q[na.index] <- NA
-        q.no.na <- q[!na.index]
-        for (i in c("p", "mean1", "sd1", "mean2", "sd2", "p.mix")) assign(i, 
-            arg.mat[!na.index, i])
-        if (any(p < 0 | p > 1)) 
-            stop("All non-missing values of 'p' must be between 0 and 1.")
-        if (any(c(sd1, sd2) < .Machine$double.eps)) 
-            stop("All non-missing values of 'sd1' and 'sd2' must be positive.")
-        if (any(p.mix < 0 | p.mix > 1)) 
-            stop("All non-missing values of 'p.mix' must be between 0 and 1.")
-        q.no.na[p == 0] <- -Inf
-        q.no.na[p == 1] <- Inf
-        index <- (1:length(q.no.na))[0 < p & p < 1]
-        if (any(index)) {
-            o.fcn <- function(q, mean1, sd1, mean2, sd2, p.mix, 
-                p) {
-                (pnormMix(q, mean1, sd1, mean2, sd2, p.mix) - 
-                  p)^2
-            }
-            for (i in index) {
-                q.no.na[i] <- nlminb(start = (1 - p.mix[i]) * 
-                  qnorm(p[i], mean1[i], sd1[i]) + p.mix[i] * 
-                  qnorm(p[i], mean2[i], sd2[i]), o.fcn, mean1 = mean1[i], 
-                  sd1 = sd1[i], mean2 = mean2[i], sd2 = sd2[i], 
-                  p.mix = p.mix[i], p = p[i])$par
-            }
-        }
-        q[!na.index] <- q.no.na
-    }
-    if (!is.null(names.p)) 
-        names(q) <- rep(names.p, length = length(q))
-    q
-}
+qnormMix <-
+function (p, mean1 = 0, sd1 = 1, mean2 = 0, sd2 = 1, p.mix = 0.5)
+{
+    if (any(p < 0 | p > 1))
+        stop("All non-missing values of 'p' must be between 0 and 1.")
+    if (any(c(sd1, sd2) < .Machine$double.eps))
+        stop("All non-missing values of 'sd1' and 'sd2' must be positive.")
+    if (any(p.mix < 0 | p.mix > 1))
+        stop("All non-missing values of 'p.mix' must be between 0 and 1.")
+    sapply(p, function(y)
+        uniroot(function(x) pnormMix(x,mean1,sd1,mean2,sd2,p.mix)-y,
+            interval = range(qnorm(y,mean1,sd1),qnorm(y,mean2,sd2)),
+            tol = 10^{-16})$root)
+}

tests/testthat/test_qnormMix.R

@@ -0,0 +1,23 @@
+library(EnvStats)
+
+
+test_that("qnormMix 1",{
+  q <- qnormMix(.95,mean1=-1,sd1=1, mean2=4,sd2=1, p.mix=.2)
+  expect_identical(round(q,2),4.67)
+})
+
+test_that("qnormMix 2",{
+  q <- qnormMix(.95,mean1=-1.5,sd1=1, mean2=6,sd2=1, p.mix=.2)
+  expect_identical(round(q,2),6.67)
+})
+
+test_that("qnormMix 3",{
+  q <- qnormMix(.95,mean1=-5,sd1=1, mean2=5,sd2=1, p.mix=.5)
+  expect_identical(round(q,2),6.28)
+})
+
+test_that("qnormMix 4",{
+  q <- qnormMix(0.5, 10, 2, 20, 2, 0.1)
+  expect_identical(round(q,2),10.28)
+})
+