Random sample from given bivariate discrete distribution

Question

Suppose I have a bivariate discrete distribution, i.e. a table of probability values P(X=i,Y=j), for i=1,...n and j=1,...m. How do I generate a random sample (X_k,Y_k), k=1,...N from such distribution? Maybe there is a ready R function like:

sample(100,prob=biprob)

where biprob is 2 dimensional matrix?

One intuitive way to sample is the following. Suppose we have a data.frame

dt=data.frame(X=x,Y=y,P=pij)

Where x and y come from

expand.grid(x=1:n,y=1:m)

and pij are the P(X=i,Y=j).

Then we get our sample (Xs,Ys) of size N, the following way:

set.seed(1000) 
Xs <- sample(dt$X,size=N,prob=dt$P)
set.seed(1000)
Ys <- sample(dt$Y,size=N,prob=dt$P)

I use set.seed() to simulate the "bivariateness". Intuitively I should get something similar to what I need. I am not sure that this is correct way though. Hence the question :)

Another way is to use Gibbs sampling, marginal distributions are easy to compute.

I tried googling, but nothing really relevant came up.

Answer 1

You are almost there. Assuming you have the data frame dt with the x, y, and pij values, just sample the rows!

dt <- expand.grid(X=1:3, Y=1:2)
dt$p <- runif(6)
dt$p <- dt$p / sum(dt$p)  # get fake probabilities
idx <- sample(1:nrow(dt), size=8, replace=TRUE, prob=dt$p)
sampled.x <- dt$X[idx]
sampled.y <- dt$Y[idx]

Answer 2

It's not clear to me why you should care that it is bivariate. The probabilities sum to one and the outcomes are discrete, so you are just sampling from a categorical distribution. The only difference is that you are indexing the observations using rows and columns rather than a single position. This is just notation.

In R, you can therefore easily sample from your distribution by reshaping your data and sampling from a categorical distribution. Sampling from a categorical can be done using rmultinom and using which to select the index, or, as Aniko suggests, using sample to sample the rows of the reshaped data. Some bookkeeping can take care of your exact case.

Here's a solution:

library(reshape)

# Reshape data to long format.
data <- matrix(data = c(.25,.5,.1,.4), nrow=2, ncol=2)
pmatrix <- melt(data)

# Sample categorical n times.
rcat <- function(n, pmatrix) {
    rows <- which(rmultinom(n,1,pmatrix$value)==1, arr.ind=TRUE)[,'row']
    indices <- pmatrix[rows, c('X1','X2')]
    colnames(indices) <- c('i','j')
    rownames(indices) <- seq(1,nrow(indices))
    return(indices)
}

rcat(3,pmatrix)

This returns 3 random draws from your matrix, reporting the i and j of the rows and columns:

  i j
1 1 1
2 2 2
3 2 2