From: Martin Maechler <maechler_at_stat.math.ethz.ch>

Date: Sat 11 Feb 2006 - 09:16:10 EST

}

r

}

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

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 veryMartin> 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 PLEASEMartin> 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
*