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

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")

• Can you clarify what "I want to generate a sample of random variates choosing 5 out of 7 days of week such that the aggregate counts resemble a given day-of-week profile" means? – dsaxton Feb 23 '16 at 21:01
• Means when counting the number of times each day of week appears in the sample, the proportions with respect to the total match the desired proportions. – arameri Feb 23 '16 at 21:13

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)
}


> 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

• Thanks Xi'an. I understand your point. The sample() function in R already adjusts the probabilities when the sampling is without replacement (according to the documentation). My custom function is trying to do the same numerically. I am not concerned about the difference of the two sampling methods. What I am concerned is about the difference between the desired counts distribution and the resulting counts distribution. – arameri Feb 23 '16 at 21:56
• Yes, the resulting imbalance is the problem. How to overcome that imbalance is the exact question still unanswered. I have modified the counts by trial-an-error to achieve approximate balance. However, I believe there must be some theoretical way to calculate. – arameri Feb 24 '16 at 22:54