# How to calculate conditional probabilities in AnyDice?

While writing the addendum to this answer, which considers the relative value of skill vs. characteristic in the "3d20 system" of Neuroshima, I found myself wanting an answer to a deceptively simple question: how many skill points are needed to succeed if the lowest roll is a natural success, vs. if it's not? In other words, I basically wanted to plot the distributions of:

• the middle roll of 3d20, given that the lowest roll is less than some given threshold x; and
• the sum of the lowest and middle rolls, given that the lowest roll is at least x.

In statistics, this would be just a standard conditional probability distribution, e.g. $$p_x(y) = P(Y = y \mid X < x),$$ $$q_x(z) = P(X + Y = z \mid X \ge x),$$ where $$\X\$$ and $$\Y\$$ are (interdependent) random variables representing the lowest and the middle roll of 3d20 respectively. You could compute this easily just by taking the joint distribution of $$\(X,Y)\$$, dropping those cases where the condition (e.g. $$\X < x\$$) fails, rescaling the remaining probabilities so that they sum to 1 and then optionally summing over the conditioning variable $$\X\$$ to obtain the marginal distribution of $$\Y\$$ (or $$\X + Y\$$).

Unfortunately, there seems to be no simple built-in way to do this in AnyDice. In fact, there doesn't even seem to be any way to answer simpler conditional probability questions like, say "what is the average sum of 3d6 if the sum rolled is even, vs. if it's odd?"

So, hence this question: Is there any way to calculate a conditional probability distribution in AnyDice, and if so, how?

Disclaimer: I realize that this question may be borderline off-topic for this site, as it's more of a programming / math question. That said, it did arise in an RPG-related context — specifically, while writing an answer here on RPG.SE — and I suspect the answer(s) may be useful to others using AnyDice to answer similar questions about other systems as well. I'll let the community decide if this Q&A should stay here or not.

Also, I did eventually manage to come up with a (slightly hacky but workable) solution to my problem on my own, so I've posted a self-answer below. That said, other answers are more than welcome too. If there's a better way to achieve this, I would very much like to know it.

It turns out that there is a way to do this in AnyDice, at least kind of. It's a bit hacky, but it works.

The secret is rerolling.

Specifically, one general way to sample from a conditional probability distribution is called rejection sampling. Basically, you sample a value from the original (unconditioned) probability distribution, and if it fails the condition, you reject it and keep resampling until you get a result that does satisfy the condition.

And we can simulate this process in AnyDice. For example, here's a simple AnyDice function that takes a die and rerolls it if its value is not in a given range:

function: restrict ROLL:n to RANGE:s else REROLL:d {
if ROLL = RANGE { result: ROLL }
else { result: REROLL }
}
function: restrict ROLL:d to RANGE:s once {
result: [restrict ROLL to RANGE else ROLL]
}


However, this only models one reroll, but that's fine. We can just iterate it:

function: restrict ROLL:n to RANGE:s else REROLL:d {
if ROLL = RANGE { result: ROLL }
else { result: REROLL }
}
function: restrict ROLL:d to RANGE:s {
loop I over {1..20} {
ROLL: [restrict ROLL to RANGE else ROLL]
}
result: ROLL
}


Now, you might look at this code and think that this still only does 20 rerolls, but that's actually not the case. Rather, it effectively does 220, or about a million rerolls! The reason for that surprising efficiency is because we update the ROLL distribution on each iteration. So on the second iteration, we're sampling from the already re-rolled distribution, and if the sample is rejected, resampling from the same already re-rolled distribution. So, basically, each iteration doubles the effective number of rerolls.

A million rerolls isn't quite infinitely many, but it's pretty close for most purposes. And if it's really not enough (which we can easily spot in the output, by the fact that the supposedly rejected values appear in it with a non-zero probability), we can always increase the iteration count from 20 to, say, 30 for a billion effective rerolls.

output [restrict 3d6 to {3,5,7,9,11,13,15,17}] named "3d6 if odd"
output [restrict 3d6 to {4,6,8,10,12,14,16,18}] named "3d6 if even"


and an example of the output:

(Surprisingly enough, it turns out that the average the same in both cases!)

But how can we use this function to handle more complex cases, like the original example of "middle of 3d20 if lowest is less than $$\x\$$", where the variable whose distribution we want is not the same as what we want to condition it on?

Well, one fairly simple way is to write a function that takes in the input roll (here, 3d20) as a sequence and maps it to the output we want (i.e. the middle roll, in this case) while also mapping any rejected cases to some bogus outcome such as −1. Then we can just use the function above to reject the bogus outcome, and get the conditional distribution we want, e.g. like this:

function: middle of ROLL:s if lowest in RANGE:s {
if 3@ROLL = RANGE { result: 2@ROLL } \ assumes a three die pool! \
else { result: -1 }
}

MAX: 10
DIST: [middle of 3d20 if lowest in {1..MAX}]

output DIST named "middle of 3d20 if lowest <= [MAX] (else -1)"
output [restrict DIST to {1..20}] named "middle of 3d20 if lowest <= [MAX] (conditional)"


The actual script I wrote for the original answer that inspired this Q&A is similar, although it makes use of an additional dice relabeling trick (described in the answer) to easily calculate the number of skill points needed to bring both of the lowest rolls under the threshold.