Dear all,
1. Let me open the question a bit: how to sample from a truncated
MULTIVARIATE normal, with known parameters. I.e., given
MVN(m,V) where m is a p-vector and V is the pxp variance matrix,
and given a p-vector a, I want to sample X=(x1,...,xp) ~ MVN(m,V)
truncated to the region {x1 < a1, x2 < a2,..., xp < ap}
The solution I've used before is a Gibbs sampler where at each step the
conditional distribution is univariate truncated normal, for which I'm
using Alan's method. Are there any viable competitors? Any special
issues one should be aware of?
2. I'll go further now and assume what I want is to compute the
PROBABILITY of a rectangular region for a MVN distribution. I can use the
simulation above for this problem, but are there better ways to handle
this, for say p >= 3?
Thanks for your insights,
Florin
On Wed, 4 Apr 2001, Alan Zaslavsky wrote:
>
> A number of fantastical solutions have been proposed to the problem of
> generating data from a truncated normal distribution. Brian Ripley
> and Nick Ellis have put forward one of the CORRECT solutions, i.e. to
> draw from an untruncated normal and reject values outside the range
> of truncation.
>
> Brian and Nick's solution is the simplest, and will be useful and even
> fast if the range of truncation is not too narrow. It begins to break
> down as you truncate to a range with small probability content under
> the untruncated distribution, since then you will have a correspondingly
> small acceptance rate.
>
> The other general approach to this problem is to generate from the
> appropriately truncated uniform and use the inverse probability function
> to transform back to the truncated range of the desired distribution.
> For example:
>
> trunc.rnorm <- function(n,mean,sd,lower,upper)
> {
> zlower <- (lower-mean)/sd ## transform bounds to z-scores
> zupper <- (upper-mean)/sd
> u <- runif(n,pnorm(zlower),pnorm(zupper)) ## truncated uniform
> mean + sd * qnorm(u) ## transform to normal
> }
>
> The computations here are slower per draw than rnorm(), but nothing is
> discarded. If you go far into the tail, the numerical approximations in
> pnorm() and qnorm() begin to break down and then you need to use other
> approximations.
>
> This method obviously generalizes to any distribution for which you have
> a "p---" and "q---" function.
> ---------------------------------------------------------------------
> This message was distributed by s-news@lists.biostat.wustl.edu. To
> unsubscribe send e-mail to s-news-request@lists.biostat.wustl.edu with
> the BODY of the message: unsubscribe s-news
>
============================
Florin Vaida
Assist Prof, Dept of Biostatistics
Harvard School of Public Health
655 Huntington Ave Boston MA 02115
Tel (617) 432-2914
|