[R] mixture normal distributions

[Martin Maechler]

>>>>> "Martin" == Martin Maechler <maechler at 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 at yahoo.com
    >>> 
    >>> Al-Eid
    >>> 
    >>> The university of Manchester.

    Martin> ______________________________________________
    Martin> R-help at 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