R help - Nov 2010 - How to generate multivariate uniform distribution random numbers?

I wish to generate 100 by 1 vector of x1 and x2 both are uniform distributed
with covariance matrix \Sigma.

Thanks,

Michael

	[[alternative HTML version deleted]]

On 06-Nov-10 21:54:41, michael wrote:
> I wish to generate 100 by 1 vector of x1 and x2 both are uniform
> distributed with covariance matrix \Sigma.
> 
> Thanks,
> Michael
First, some comments.

1. I don't think you mean a "100 by 1 vector of x1 and x2" since
   you have two variables. "100 by 2" perhaps?

2. Given the range over which x1 is to be uniformly distributed
   (say of length A) and the range for x2 (say of length B),
   then by virtue of the uniform distribution the variance of x1
   will be (A^2)/12 and the variance of x2 will be (B^2)/12,
   so you don't have a choice about these. So your only "free
   parameter" is the covariance (or the correlation) between
   x1 and x2.

That said, it is a curious little problem. Here is one simple
solution (though it may have properties that you do not want).
Using X and Y for x1 and x2 (for simplicity), and using ranges
(-1/2,1/2) for the ranges of X and Y (you could rescale later):

1. Generate X uniformly distributed on (-1/2,1/2).
2. With probability p (0 < p < 1) let Y=X.
   Otherwise, let Y be independently uniform on (-1/2,1/2).

This can be implemented by the following code:

  n  <- 100
  p  <- 0.6            # (say)
  X  <- runif(n)-1/2
  Z  <- runif(n)-1/2
  U  <- rbinom(n,1,p)  # =1 where Y=X, = 0 when Y independent of X
  Y  <- X*U + Z*(1-U)
  XY <- cbind(X,Y)     # your n by 2 matrix of bivariate (X,Y)

Now the marginal distributions of X and Y are both uniform.
var(X) = 1/12 and var(Y) = 1/12. Since the means are 0,
cov(X,Y) = Exp(X*Y) = p*Exp(X^2) = p/12
cor(X,Y) = cov(X,Y)/sqrt(var(X)*var(Y)) = (p/12)/(1/12) = p.

So all you need to do to get a desired correlation rho between
X and Y is to set p = rho.

The one thing you may not like about this is the diagonal line
you will get for the cases where X=Y:

  plot(X,Y)

Test:

  n <- 100000
  p  <- 0.6 # (say)
  X  <- runif(n)-1/2
  Z  <- runif(n)-1/2
  U  <- rbinom(n,1,p)
  Y  <- X*U + Z*(1-U)
  var(X)
  # [1] 0.08340525  # theory: var = 1/12 = 0.08333333
>   var(Y)
[1] 0.08318914      # theory: var = 1/12 = 0.08333333
>   cov(X,Y)
[1] 0.04953733      # theory: cov = p/12 = 0.6/12 = 0.05
>   cor(X,Y)
[1] 0.5947063       # theory: cor = p    = 0.6

It would be interesting to see a solution which did not involve
having cases with X=Y!

Hoping this helps,
Ted.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <ted.harding at wlandres.net>
Fax-to-email: +44 (0)870 094 0861
Date: 06-Nov-10                                       Time: 23:07:26
------------------------------ XFMail ------------------------------