Generating n distinct combinations from a weighted distribution If I'm not mistaken, this is an extension of this question:
How to resample in R without repeating permutations?
I have a vector x representing a non-uniform discrete distribution (with probability vector p where sum(p)==1) from which I want to generate $m$ distinct draws of size $k$.
For $\frac{m}{\binom{ || x || }{k}}$ very small, a simple while loop basically suffices:
# have to sort or unique() won't equate (3,4) and (4,3)
do_sample = function(mm)
  replicate(mm, sort(sample(x, k, prob = p)), simplify = FALSE)
draws = unique(do_sample(m))
while (length(draws) < m) {
  draws = unique(c(draws, do_sample(m - length(draws))))
}

IINM this is a more readable but inefficient version of @whuber's suggestion -- calculating unique all the time gets expensive quickly.
For $\frac{m}{\binom{ || x || }{k}}$ large (e.g. $||x||=72$, $m=17000$, $k=4$, this becomes nigh interminable as the probability of drawing a distinct combination appears to plummet (the distribution is somewhat skewed, so the "rare" draws are increasingly difficult to happen upon)
Given that $\binom{72}{4}$ is manageably small, in the uniform case I would just do 
universe = combn(x, 4, simplify = FALSE)
universe[sample(length(universe), m)]

But the different combinations don't have equal probabilities on account of p being non-uniform, so this doesn't work.
The closed-form probability of drawing $(x_1, x_2, x_3)$ is (IINM):
$$ \tilde{\pi}\left( x_1, x_2, x_3 \right) = \pi(x_1) \pi(x_2) \pi(x_3) \left( \frac{1}{1-\pi(x_1)} \frac{1}{1 - \pi(x_1) - \pi(x_3)} + \frac{1}{1-\pi(x_1)} \frac{1}{1 - \pi(x_1) - \pi(x_2)} + \frac{1}{1-\pi(x_2)} \frac{1}{1 - \pi(x_2) - \pi(x_3)} + \frac{1}{1-\pi(x_2)} \frac{1}{1 - \pi(x_2) - \pi(x_1)} + \frac{1}{1-\pi(x_3)} \frac{1}{1 - \pi(x_3) - \pi(x_1)} + \frac{1}{1-\pi(x_3)} \frac{1}{1 - \pi(x_3) - \pi(x_2)}\right) $$
A similar version would exist for $\tilde{\pi}\left( x_1, x_2, x_3, x_4 \right)$, and these could be used as the prob argument in sample(length(universe), m, prob=), but these are pretty tedious to write up; I can't help but wonder if there's a better way to go about this.

Editing to add a reasonably concise way to calculate the probabilities:
# permutations gets all permutations of 1:n, see e.g.
#   https://stackoverflow.com/a/20199902/3576984
# the sum on the RHS of \tilde{\pi} has k! terms,
#   one corresponding to each (ordered) arrangement of
#   the k chosen objects
all_perm = permutations(k)
prob = function(k) {
  perm = permutations(k)
  combn(seq_along(x), k, function(idx) {
    p_idx = p[idx]
    # an array with the n! arrangements of the k probabilities;
    #   the last column is not needed, so drop it
    p_perm = array(p[perm[ , -k]], dim = dim(perm)-0:1)
    # e.g., for k=3, this is 
    #   1-rowSums(p_perm[ , 1])
    #   1-rowSums(p_perm[ , 1:2])
    # similar to rowCumSums (but apply requires iterating over rows-->slow)
    prod(p_idx) * sum(1/Reduce(`*`, lapply(1:(k-1L), function(J)
      1 - rowSums(p_perm[ , 1:J, drop = FALSE]))))
  })
}

 A: Could you say a bit more about how you derived the joint distribution? 
If I understand correctly, prob should output a valid PMF, but it doesn't seem to add to unity:
# Original suggested solution
permutations <- function(n){
  if(n==1){
    return(matrix(1))
  } else {
    sp <- permutations(n-1)
    p <- nrow(sp)
    A <- matrix(nrow=n*p,ncol=n)
    for(i in 1:n){
      A[(i-1)*p+1:p,] <- cbind(i,sp+(sp>=i))
    }
    A
  }
}

prob <- function(k, x, p) {
  perm <- permutations(k)
  combn(seq_along(x), k, function(idx) {
    p_idx <- p[idx]
    p_perm <- array(p[perm[, -k]], dim = dim(perm) - 0:1)
    prod(p_idx) * sum(1 / Reduce(`*`, lapply(1:(k-1L), function(J)
      1 - rowSums(p_perm[, 1:J, drop = FALSE]))))
  })
}

# Parameters
lambda <- 5
x <- 1:12
k <- 3
p <- dpois(x, lambda) / sum(dpois(x, lambda))

combinations <- as.data.frame(t(combn(x, k)))
probs <- prob(k, x, p)
sum(probs)  # 0.8649801

# Alternative solution
get_prob <- function(s, p) {

  get_prob_for_perm <- function(s, p) {
    if(length(s) <= 0) return(1)
    x <- s[1]
    p_x <- p[x]
    p[x] <- 0
    p_x * get_prob_for_perm(s[-1], p / sum(p))
  }

  perm <- combinat::permn(s)
  sum(sapply(perm, function(s) get_prob_for_perm(s, p)))
}

combinations$p <- apply(combinations, 1, function(s) get_prob(s, p))
sum(combinations$p) # 1

# Sampling
combinations[sample(nrow(combinations), k, prob = probs2), ]

I'm calculating the probabilities quite naively above, so I'm not sure they're correct.
