title: Random sample from given bivariate discrete distribution
body: <p>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:</p>
  
  <pre><code>sample(100,prob=biprob)
  </code></pre>
  
  <p>where biprob is 2 dimensional matrix? </p>
  
  <p>One intuitive way to sample is the following. Suppose we have a data.frame</p>
  
  <pre><code>dt=data.frame(X=x,Y=y,P=pij)
  </code></pre>
  
  <p>Where x and y come from</p>
  
  <pre><code>expand.grid(x=1:n,y=1:m)
  </code></pre>
  
  <p>and pij are the P(X=i,Y=j). </p>
  
  <p>Then we get our sample (Xs,Ys) of size N, the following way:</p>
  
  <pre><code>set.seed(1000) 
  Xs &lt;- sample(dt$X,size=N,prob=dt$P)
  set.seed(1000)
  Ys &lt;- sample(dt$Y,size=N,prob=dt$P)
  </code></pre>
  
  <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 :) </p>
  
  <p>Another way is to use Gibbs sampling, marginal distributions are easy to compute. </p>
  
  <p>I tried googling, but nothing really relevant came up.</p>