[R] simulate from conditional distribution

[Petr Savicky]

On Thu, Sep 13, 2012 at 05:02:54PM -0400, li li wrote:
> Dear all,
>         Y, X are bivariate normal with all the parameters known.
> I would like to generate numbers from the distribution Y | X > c
> where c is a constant.
>         Does there exist an R function generating
> random numbers from such a distribution?

Hi.

One of the options, how to generate such numbers, is the rejection
method. This works well, if Pr(X > c) is not too small, but may be
very bad otherwise.

For generating a single number, try the following

  library(mvtnorm)
  sigma <- rbind(c(3, 1), c(1, 2))
  cc <- 2
  while (1) {
      v <- rmvnorm(1, sigma=sigma)
      if (v[2] > cc) break
  }
  x <- v[1]

For a larger number, try something like the following

  n <- 500
  for (k in 2^(1:10)) {
      v <- rmvnorm(k*n, sigma=sigma)
      x <- v[v[, 2] > cc, 1]
      if (length(x) >= n) break
  }
  stopifnot(length(x) >= n)
  x <- x[1:n]

Hope this helps.

Petr Savicky.