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.
>>>
>>> Many thanks for your help in advance.
>>>
>>> My E-mail is aleid2001@yahoo.com
>>>
>>> Al-Eid
>>>
>>> The university of Manchester.

    Martin> ______________________________________________
    Martin> R-help@stat.math.ethz.ch mailing list
    Martin> https://stat.ethz.ch/mailman/listinfo/r-help PLEASE
    Martin> do read the posting guide!
    Martin> http://www.R-project.org/posting-guide.html

R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Sat Feb 11 09:19:35 2006

This archive was generated by hypermail 2.1.8 : Mon 13 Feb 2006 - 22:58:43 EST