How to generate a more accurate distribution of sample of random variates?

Question

I want to generate a sample of random variates choosing 5 out of 7 days of the week, $X_1,\ldots,X_5$, such that the aggregate counts resemble a given day-of-week profile, namely $$\mathbb{E}[\mathbb{I}_i(X_1)+\cdots+\mathbb{I}_i(X_5)]=5p_i$$ for $i=1,\ldots,7$.

I am using the traditional inverse transform method with pseudo-random uniform U numbers and the desired cumulative distribution. Why is the sample distribution flatter than the generating distribution?

I am assuming it has to do somehow with conditional probabilities. Assuming the pseudo-random numbers are acceptable, what is the best way to accomplish this?

Below stands my implementation in R:

counts <- c(2268, 3232, 4173, 4113, 4079, 3877, 2648)
dow.pdf <- counts / sum(counts)
dow.cdf <- cumsum(dow.pdf)

randDow <- function(n) {
  dow <- NULL
  while (length(dow) < n) {
    #generate pseudo-random number
    U <- runif(1, min = 0, max = 1)

    #find the position in the CDF
    rv <- which.min(U > dow.cdf)

    #check for duplicates
    if (!(rv %in% dow)) {
      dow <- c(dow, rv)
    }
  }
  return(sort(dow))
}

#generate 10000 random variates method 1
s1 <- t(replicate(10000,c(sort(sample(1:7,5,rep=F,prob=dow.pdf)))))

#generate 10000 random variates method 2
s2 <- t(replicate(10000, randDow(5)))

#summarize sample results
ps1 <- table(s1)/sum(table(s1))
ps2 <- table(s2)/sum(table(s2))

#compare visually
barplot(matrix(c(dow.pdf, ps1, ps2), nr=3, byrow=TRUE), beside=T, 
    names.arg=1:7, main="Day of Week Distribution")
legend("topleft", c("pdf","ps1", "ps2"), pch=15, cex=0.7,
   bty="n")

Answer 1

Update [02/30]: As my colleague Robin Ryder pointed out to me, the problem is solved in Brewer's book Sampling with unequal probabilities. This resolution is actually discussed in an earlier Cross validated question. It is however unclear to me whether or not there exists a sampling strategy with fixed probabilities, in the sense of running

sample(1:N,n,rep=FALSE,prob=q)

to produce the sample with expected frequencies $p_i$...

---

The traditional inverse cdf method works for a single simulation. However, if you generate 5 draws without replication, according to the frequencies given by counts, the marginal distribution of those draws is no longer the empirical cdf dow.cdf. The fact that both your methods lead to highly similar frequencies is reflecting upon this difference.

For instance, consider the formal case when counts=c(1000,100,1,1)) and you draw n=2 out of 4 without replacement: because of the huge imbalance between the first two counts and the last two counts, most draws by sample((1:4),2,rep=FALSE,prob=c(1e3,1e2,1,1)) will be 1 2 despite 2 having a probability ten times smaller. As shown by the picture below:

However, if the number of draws is reasonably small against the number of possible values and if the probabilities are not too different, there may exist a solution. For instance, if one generates only 2 draws out of 7 possible values, and if one wants to achieve the probability vector $\mathbf{q}$ on the resulting sample, the sampling probability vector $\mathbf{p}$ is solution of the equations $(i=1,\ldots,7)$ $$ 2q_i=p_i+\sum_{j\ne i} p_j\,\frac{p_i}{1-p_j} $$ While solving this equation in a closed form is not feasible, a recursion based on the fact that the uniform probability is a fixed point may produce a solution. Here is an example where it works:

#function to be inverted
targ=function(p){
  p*(2-length(p)+sum(1/(1-p))-1/(1-p))}
#component in the above
nores=function(p,i=1){
  2-length(p)+sum(1/(1-p[-i]))}
#target probability
kuh=(1:7)/sum((1:7))
#probability when resampling 
peh=rep(1,7)/7
while (sum(abs(2*kuh-targ(peh)))>.001){
  i=sample(1:7,1)
  peh[i]=2*kuh[i]/nores(peh,i)
  peh=peh/sum(peh)
  }

Running this code leads to

> peh
[1] 0.03261515 0.06597545 0.10094319 0.13761820 0.17673115 0.21921258 0.26690427

and checking that it works:

> s1 <- t(replicate(10000,c(sort(sample(1:7,2,rep=F,prob=peh)))))
> table(s1)/sum(table(s1))
s1
      1       2       3       4       5       6       7 
0.03710 0.07225 0.10695 0.13990 0.17630 0.21545 0.25205 
> kuh
[1] 0.03571429 0.07142857 0.10714286 0.14285714 0.17857143 0.21428571 0.25000000

But if one picks a more extreme repartition of the draws, like

kuh=(1:7)^5/sum((1:7)^5)

there is no solution (found by this recursive approach). Extending to more draws is equally feasible if a wee bit more cumbersome. For instance, here is the version for 3 draws:

targ=function(p){
  q=rat=p/(1-p)
  for (i in 1:7){
    dble=outer(rat[-i],p[-i],'*')
    tple=dble/(1-outer(p[-i],p[-i],'+'))
    q[i]=1+sum(rat[-i])+sum(tple)-sum(diag(tple))}
  return(p*q)}
nores=function(p,i=1){
  rat=p/(1-p)
  dble=outer(rat[-i],p[-i],'*')
  tple=dble/(1-outer(p[-i],p[-i],'+'))
  1+sum(rat[-i])+sum(tple)-sum(diag(tple))}

which leads to a solution using the iterative pluggin approach as follows:

kuh=(1:7)/sum((1:7))
peh=rep(1,7)/7
while (max(abs(3*kuh-targ(peh)))>.001){
  i=sample(1:7,1)
  peh[i]=3*kuh[i]/nores(peh,i)
  peh=peh/sum(peh)}

as it stops at

> peh
[1] 0.02797048 0.05788001 0.09031821 0.12651015 0.16894929 0.22434742 0.30402444

which is close enough to the target:

> max(abs(3*kuh-targ(peh)))
[1] 0.0004822188