# Re: [R] mixture normal distributions

From: Martin Maechler <maechler_at_stat.math.ethz.ch>
Date: Sat 11 Feb 2006 - 09:16:10 EST

>>>>> "Martin" == Martin Maechler <maechler@stat.math.ethz.ch> >>>>> on Fri, 10 Feb 2006 22:20:11 +0100 writes:

>>> Dear R helper, I mange to transform uniform sequences to
>>> mixture normal distributions using the following cods:

```    Martin> you forgot to mention the important fact that you
Martin> are working with package "nor1mix" (of which I am
Martin> the maintainer which you could have seen from
Martin> library(help = nor1mix) or
Martin> packageDescription("nor1mix")

```

>>> > K<-50000 > prime<-c(29) , where 29 is prim number >
>>> UN<-seq(1:K)%*%t(sqrt(prime)) > U1<-UN-as.integer(UN) >
>>> e<-norMix(mu=c(-0.825,0.275), sig2 = c(0.773,0.773), w =
>>> c(0.25,0.75), name = NULL, long.name = FALSE) >
>>> U<-matrix(qnorMix(e,U1),K,1),

Martin> This looks like you want to generate pseudo-random     Martin> values from the mixture distribution.

```    Martin> However by the slowest most possible way: qnorMix()
Martin> is really slow because it calls uniroot() for each
Martin> value.

Martin> rnorMix() will generate random values distributed
Martin> according to the normal mixture very very
```
Martin> efficiently (particularly compared to using     Martin> qnorMix()).
```    Martin> But indeed, you have detected a bug in qnorMix()
Martin> (that happens pretty rarely).  Indeed your example
Martin> also shows that the algorithm can be much improved
Martin> in some cases.

```

Martin> The next verion of nor1mix should contain a more     Martin> "robust" qnorMix() function.

in the mean time, you can use the following which already fixes the problem you found :

qnorMix <- function(obj, p)
{
if (!is.norMix(obj))

stop("'obj' must be a 'Normal Mixture' object!")   mu <- obj[, "mu"]
sd <- sqrt(obj[, "sig2"])
k <- nrow(obj)# = #{components}
if(k == 1) # one component

return(qnorm(p, mu, sd))

## else

## vectorize in `p' :
r <- p
left <- p <= 0 ; r[left] <- -Inf
right <- p >= 1 ; r[right] <- Inf
imid <- which(mid <- !left & !right) # 0 < p < 1   ## p[] increasing for easier root finding start:   p <- sort(p[mid], index.return = TRUE)   ip <- imid[p\$ix]
pp <- p\$x
for(i in seq(along=pp)) {

```      rq <- range(qnorm(pp[i], mu, sd))
## since pp[] is increasing, we can start from last 'root':
if(i > 1) rq[1] <- root
## make sure, 'lower' is such that f(lower) < 0 :
delta.r <- 0.01*abs(rq[1])
ff <- function(l) pnorMix(obj,l) - pp[i]
while(ff(rq[1]) > 0) rq[1] <- rq[1] - delta.r

root <- uniroot(ff, interval = rq)\$root
r[ip[i]] <- root
```

}
r
}

I hope this helps you.

>>> But somtimes if i use ,e.g, 23 or 11 instead of 29 it
>>> will give me the following error.
>>>
>>> > K<-30000 > prime<-c(23) >
>>> UN<-seq(1:K)%*%t(sqrt(prime)) > U1<-UN-as.integer(UN) >
>>> e<-norMix(mu=c(-0.825,0.275), sig2 = c(0.773,0.773), w =
>>> c(0.25,0.75), name = NULL, long.name = FALSE) >
>>> U<-matrix(qnorMix(e,U1),K,1) Error in
>>> uniroot(function(l) pnorMix(obj, l) - pp[i], interval =
>>> rq) : f() values at end points not of opposite sign
>>>
>>>
>>> I am seeking help how to avoid this error.
>>>
>>>
>>> My E-mail is aleid2001@yahoo.com
>>>
>>> Al-Eid
>>>
>>> The university of Manchester.

```    Martin> ______________________________________________
Martin> R-help@stat.math.ethz.ch mailing list