If I understand the question, and I may not, you already know the n different pivalues, you want to draw new samples where you learn the result of each of the n experiments, but you are only interested in those sets of results which have nr "red" and n-nr"blue".  A lazy way is to draw a set of samples and check whether the set has exactly nr "red". If not, discard the results and repeat. This is efficient in terms of your time spent programming, but probably not in terms of computer resources.

Or, when nr < n, distribute nrreds among the n experiments, keeping the relative probabilities of red. Pick one experiment to get a red by drawing a sample of size 1 from a categorical probability distribution with probabilities Pi = pi / (sum(pi) ). Change the corresponding pi to 0, recalculate Pivalues, and repeat for a total of nrdraws.

I just want to add few words to suggestions given by Charles. 1) To "A lazy way". Use binary representations for numbers from 0 to 2^n-1 and attach the relative probabilities to those numbers whose binary representations have exactly n_r (or n-n_r) ones (and 0 to the remaining numbers) . For instance, when n=4 and n_r=2 non-trivial probabilities have the numbers 3=0011(bin), 5=0101(bin), 6=0110(bin), 9=1001(bin), 10=1010(bin), 12=1100(bin).  Then use random variable X with random values ranging from 1 to {n\choose n_r} and  calculated relative probabilities. In the above example, P(X=1)=P(result=0011), ...,  P(X=6)=P(result=1100). 

2) Or, you can take 2 counters, say, R(red) and B(blue). At the beginning R=B=0. Then each time t (1

Dear Ton,

I observe that you really want to draw nr reds AND  (n - nr) blues.

This can become a hell of a problem in general. For instance, if all the probabilities p(i) lie between 0.95 and 1,0 and you want a total of  n = 100 draws, ONE of which is red, you may be obliged to draw a Hubble time and not succeed. In order to find a viable method at all, along the lines of Yuri or Charles, the nr reds must be pretty much be the number you expect

Obviously that number is E = SUM p(k);

the variance on that number is V = SUM p(k) q(k)

So in general I expect severe, practically insurmountable problems if the nr desired is many standard deviations away from the expected.

Maybe we can talk about it in more detail tonight at the chess club! :-)

The problem is subtle, but not as severe or practically insurmountable as Manuel suggests.

My first suggestion above would draw from all possible outcomes. I noted this may not be efficient in terms of computer resources.  Manuel is correct that there are some situations where the computer may run a very long time before finding even one solution by this method.

The alternative solution I proposed will give a correct solution in nr draws even for Manuel's case (100 experiments, each pi > 0.95, final result nr = 1 and nblue=99), but it takes a bit more programming. Start with the values of pi = [0.95, 0.96, 0.99, ..., 0.95). Find the sum [97.0, for example). Calculate P values (P is not the same as p) by division, Pi = [ 0.009793814, 0.009896907, 0.01020619, ..., 0.009793814]. Draw one red and record which 'experiment' or category it was. For example, say it was i=1. We have already solved Manuel's case.  If the number of red was greater than one, you would continue. For i=1 change pi to 0, so now our working pi = [0.00, 0.96, 0.99, ..., 0.95). Find the sum [now 96.05]. Calculate new P, then continue this cycle until you have identified the nr 'experiments' that came up red. All the other 'experiments' are 'blue' since each 'experiment' results in either red or blue.  In short, you condition on or consider only those possible outcomes that end up the required number of red.

To understand the approach without the math--In the first draw, each of the 100 'experiments' has a fair chance of being chosen red. At the next draw, each of the remaining 'experiments' has a fair chance of being red, given that one 'red' has already been taken and one 'experiment' is out of the picture (p1=0 in the example).

Yuri's first method also conditions on and considers only the possible outcomes that end up with the required number of red. It begins by identifying all 2^n - 1 possible outcomes, but he then (conditions) and reduces this to the n choose nr outcomes with nr red. Then he scales the probabilities of each remaining allowable outcomes so they sum to 1.0 and randomly picks the outcome. One might run into computer memory issues at either of the first two steps if 2^n - 1 or n choose nr are large. 

Yuri's second method of conditioning may be great--I just don't understand it yet.

My second suggestion consisted in reducing  each time the pair (n_current,n_r_current)  either to the pair (n_current-1,n_r_current)  if the result of the current draw is 0, i.e. "blue", or to the pair (n_current-1,n_r_current-1)  provided the result of the current draw is 1, i.e. "red", until one gets either (n',0): "no more reds", or (n',n'): "no more blues". Just a kind of recursion.

I added also the sketch  how to recalculate the probabilities each time for the current draw.

Thank you for your answers. The "lazy approach" won't work in my case as my n is large (in some cases almost 17 million). I suppose recalculating the pi works (either in the way Charles or Yuri proposed).

My actual problem is in fact more complicated: I have more than two colors. Perhaps I'll ask that question at a later moment,

The problem arose from a practical situation. I have a dataset with data on persons in which one categorical variable is missing and I want fill in estimates for that variable for all persons (imputation). At the same time I do know the total number per category from another data source. This explains my 17 million above: there are almost 17 million inhabitants in the Netherlands.

The following method is not exact, but it is efficient and will give a very good approximation in particular when n is large. 

According to a version of the conditional limit theorem the variables xi will be approximately independent given their sum and the distribution of xi will be the element in the exponential family that gives the observed sum as its mean. That means that one should determine a parameter β such that when each of the distributions  of xi are exponentially tilted (multiplied by exp(β xi) and normalized) then their means sum up to the observed value. The distribution of xi given the sum is approximately equal to the mean after xihas been exponentially tilted. According to the law of large numbers sampling independently from these exponentially tilted distributions will give a sum that is very near the observed value. 

If one want the sampled sum to be idential with the observed sum one can modify the procedure slightly. First x1 is sampled according to the exponentially tilted distribution. Then the sample size is decreased by one and the observed number of blue balls is corrected according to the value of x1. Then x2 is sampled as before with the corrected values of the samplesize and observed number of blue balls. Etc.

This procedure will work equally well if there are more that two colors.

Dear Peter,

Thank you for your answer. I'm not sure, though, why your exact approach works if there are more than two colors. If you draw each xi separately, it's not trivial to check if you can still draw the appropriate number of balls of each color. It leads to the so-called Harem problem (a generalization of the Marriage problem known from Combinatorics).

You are right that one has a harem problem. I did not think of it as a harem problem, so obviously something appear to be missing in my argument. After thinking about it the situation appear to be better than indicated in my first answer.

First calculate the exponentially tilted distribution such that the mean values with respect to the tilted distribution equals the observed number of balle in each color. This gives a probability vector for the first ball, for the second ball etc. These probability vectors can be written columns in a matrix. Thus the column sums are 1 and the row sums are the observed number of balls of each color.

Next, consider the convex set of all matrices with entries in [0,1] such that the support of each column is a subset of the support of the original probability vectors and such that the row sums and the column sums equal the similar sums for the original matrix. It is a nice exercise to show that the extreme points of this convex set is a matrix with only zeros and ones. Therefore every color that has a positive probability for the first ball can be assigned without destroying the harem problem for the rest of the balls. 

Dear Peter,

Thank you again for your answer. I'll have to think about it some more before I fully understand it (and can prove the nice exercise).