The "uniform" in uniform prior doesn't just mean that hits and misses are equally likely. It means that you assume that you have a probability measure on the rates $[0,1]$ and this measure is the uniform measure. For example, it means the chance the true rate is between $0.9$ and $1.0$ is $0.1$. There are other measures on $[0,1]$ which have mean value $1/2$.
For example, if you start with the uniform prior and then observe $1$ hit and $1$ miss, the updated distribution is more concentrated around a rate of $1/2$ than the uniform prior. Instead of a probability of $\Delta r $ for the rate to be between $r$ and $r+\Delta r$, you would estimate the probability to be about $6 r (1-r) \Delta r$. (The factor of $6$ makes the total measure $1$.) This new distribution would also predict that the probability of getting a hit next time is $1/2$. If you observe an additional $h$ hits and $m$ misses, the expected probability of a hit is
$$\frac{h+2}{h+m+4}.$$
This is closer to $1/2$ than $\frac{h+1}{h+m+2}$. This agrees with the fact that you started with a distribution which was more concentrated near $1/2$. Of course, if you include that first hit and miss then there have been $h+1$ hits and $m+1$ misses in total.
The distribution you get from the uniform prior by observing some number of hits and misses is called a beta distribution. The uniform distribution on $[0,1]$ is $\text{Beta}(1,1)$. Beta distributions have the nice property that if you update a beta distribution with an observation, the result is still in that family. From one observation, the $\text{Beta}(a,b)$ distribution is updated either to $\text{Beta}(a+1,b)$ or $\text{Beta}(a,b+1)$. The mean of a $\text{Beta}(a,b)$ distribution is $\frac{a}{a+b}$. It is perfectly reasonable to have a prior distribution which is not a beta distribution, but the answers might not come out as nicely.
If I'm very confident that the hit rates are 50/50, I could add 2 instead of 1
to the hits and misses. Or if I'm less confident, I could add 1/2.
That might or might not agree with confidence. That is about whether you think it is likely that the rate is close to $1/2$. In some situations, you may believe the rate should be close to $0$ or close to $1$, and the first trial should be very informative. You might believe that half of your students know how to solve a type of problem, and that if you give a random student $5$ of these problems, the student is very likely to solve all $5$ or fail to solve any. This might be approximated by a $\text{Beta}(\epsilon,\epsilon)$ distribution so that after $h$ correct and $m$ incorrect, the average value is
$$\frac{h+\epsilon}{h+m+2\epsilon},$$
which is close to $0$ or $1$ after one observation.
What you really want to know is how to do this calculation quickly--preferably in your head--during game play. To that end, consider using Normal approximations to the distributions.
Using Normal approximations, we can easily determine that two rolls for damage with a 2d6+2 have about a $30\%$ chance of equalling or exceeding $20$ and three rolls for damage have about a $95\%$ chance.
Using a Binomial distribution, we can estimate there is about a $108/343$ chance of rolling twice for damage and $27/343$ chance of rolling three times. Therefore, the net chance of equaling or exceeding $20$ is approximately
$$(0.30 \times 108 + 0.95 \times 27) / 343 \approx (32 + 25)/343 = 57/343 \approx 17\%.$$
(Careful consideration of errors of approximation suggested, when I first carried this out and did not know the answer exactly, that this number was likely within $2\%$ of the correct value. In fact it is astonishingly close to the exact answer, which is around $16.9997\%$, but that closeness is purely accidental.)
These calculations are relatively simple and easily carried out in a short amount of time. This approach really comes to the fore when you just want to know whether the chances of something are small, medium, or large, because then you can make approximations that greatly simplify the arithmetic.
Details
Normal approximations come to the fore when many activities are independently conducted and their results are added up--exactly as in this situation. Because the restriction to nonnegative health (which is not any kind of a summation operation) is a nuisance, ignore it and compute the chance that the opponent's health will decline to zero or less.
There will be three rolls of the 1d20 and, contingent upon how many of them exceed the opponent's armor, from zero to three rolls of the 2d6+2. This calls for two sets of calculations.
Approximating the damage distribution. We need to know two things: its mean and variance. An elementary calculation, easily memorized, shows that the mean of a d6 is $7/2$ and its variance is $35/12 \approx 3$. (I would use the value of $3$ for crude approximations.) Thus the mean of a 2d6 is $2\times 7/2 = 7$ and its variance is $2\times 35/12 = 35/6$. The mean of a 2d6+2 is increased to $7+2=9$ without changing its variance.
Therefore,
- One roll for damages has a mean of $9$ and a variance of $35/6$. Because the largest possible damage is $14$, this will not reduce a health of $20$ to $0$.
- Two rolls for damages have a mean of $2\times 9=18$ and a variance of $2\times 35/6=35/3\approx 12$. The square root of this variance must be around $3.5$ or so, indicating the health is approximately $(20-18)/3.5\approx 0.6$ standard deviations above the mean. I might use $0.5=1/2$ for a crude approximation.
- Three rolls for damages have a mean of $27$ and variance of $35/2\approx 18$ whose square root is a little larger than $4$. Thus the health is around $1.5$ to $2$ standard deviations lower than the mean.
The 68-95-99.7 rule says that about $68\%$ of the results lie within one SD of the mean, $95\%$ within two SDs, and $99.7\%$ within three SDs. This information (which everyone memorizes) is on top of the obvious fact that no results are less than zero SDs from the mean. It applies beautifully to sums of dice.
Crudely interpolating, we may estimate that somewhere around $40\%$ or so will be within $0.6$ SDs of the mean and therefore the remaing $60\%$ are further than $0.6$ SDs from the mean. Half of those--about $30\%$--will be below the mean (and the other half above). Thus, we estimate that two rolls for damage has about a $30\%$ chance of destroying the enemy.
Similarly, it should be clear that when the mean damage is between $1.5$ and $2$ standard deviations above the health, destruction is almost certain. The 68-95-99.7 rule suggests that chance is around $95\%$.
This figure plots the true cumulative distributions of the final health (in black), their Normal approximations (in red), and the true chances of reducing the health to zero or less (as horizontal blue lines). These lines are at $0\%$, $33.6\%$, and $96.4\%$, respectively. As expected, the Normal approximations are excellent and so our approximately calculated chances are pretty accurate.
Estimating the number of rolls for damages. The comparison of a 1d20 to the armor class has three outcomes: doing nothing with a chance of $11/20$, rolling for half damages with a chance of $1/20$, and rolling for full damages with a chance of $8/20$. Tracking three outcomes over three rolls is too complicated: there will be $3\times 3\times 3=27$ possibilities falling into $10$ distinct categories. Instead of halving the damages upon equalling the armor, let's just flip a coin then to determine whether there will be full or no damages. That reduces the outcomes to an $11/20 + (1/2)\times 1/20 = 23/40$ chance of doing nothing and a $40/40 - 23/40 = 17/40$ chance of rolling for damages.
Since this is intended to be done mentally, note that the $23/40 = 8/20 + (1/2)\times 1/20 = 0.425$ is easily calculated and this is extremely close to a simple fraction $3/7 = 0.42857\ldots.$ We have placed ourselves in a situation equivalent to rolling an unfair coin with $3/7$ chance of success. This has a Binomial distribution:
- We can roll for damages twice with a chance of $3\times (4/7)\times (3/7)^2= 108/343.$
- We will roll for damages three times with a chance of $(3/7)^3 = 27/343.$
(These calculations are very easily learned; all introductory statistics courses cover the theory and offer lots of practice with them.)
Code
To verify this result (which was obtained before many of the other answers appeared), I wrote some R code to carry out such calculations in very general ways. Because they can involve nonlinear operations, such as comparisons and truncation, they do not capitalize on the efficiency of convolutions, but just do the work with brute force (using outer products). The efficiency is more than adequate for smallish distributions (having only a few hundred possible outcomes, more or less). I found it more important for the code to be expressive so that we, its users, could have some confidence that it correctly carries out what we want. Here for your consideration is the full set of calculations to solve this (somewhat complex) problem:
round <- conditional(sign(hit-armor), list(nothing, half(damage), damage))
x <- health - rep(round, n.rounds) # The battle
x <= nothing # Outcome distribution
The output is
FALSE TRUE
0.8300265 0.1699735
showing a 16.99735% chance of success (and 83.00265% chance of failure).
Of course, the data for this question had to be specified beforehand:
hit <- d(1, 20, 4) # Distribution of hit points
damage <- d(2, 6, 1) # Distribution of damage points
n.rounds <- 3 # Number of attacks
health <- as.die(20) # Opponent's health
armor <- as.die(16) # Opponent's armor
nothing <- as.die(0) # Result of no hit
This code reveals that the calculations are lurking in a class I have named die. This class maintains information about outcomes ("value") and their chances ("prob"). The class needs some basic support for creating dice and displaying their values:
as.die <- function(value, prob) {
if(missing(prob)) x <- list(value=value, prob=1)
else x <- list(value=value, prob=prob)
class(x) <- "die"
return(x)
}
print.die <- function(d, ...) {
names(d$prob) <- d$value
print(d$prob, ...)
}
plot.die <- function(d, ...) {
i <- order(d$value)
plot(d$value[i], cumsum(d$prob[i]), ylim=c(0,1), ylab="Probability", ...)
}
rep.die <- function(d, n) {
x <- d
while(n > 1) {n <- n-1; x <- d + x}
return(x)
}
die.normalize <- function(value, prob) {
i <- prob > 0
p <- aggregate(prob[i], by=list(value[i]), FUN=sum)
as.die(p[[1]], p[[2]])
}
die.uniform <- function(faces, offset=0)
as.die(1:faces + offset, rep(1/faces, faces))
d <- function(n=2, k, ...) rep(die.uniform(k, ...), n)
This is straightforward stuff, quickly written. The only subtlety is die.normalize, which adds the probabilities associated with values appearing more than once in the data structure, keeping the encoding as economical as possible.
The last function is noteworthy: d(n,k,a) represents the sum of n independent dice with values $1+a, 2+a, \ldots, k+a$. For instance, a 2d6+2 can be considered the sum of two d6+1 distributions and is created via the call d(2,6,1).
The heart of the code is the overloading of arithmetic operations. I implemented only those needed for this calculation, but did so in a way that is easy to extend, as should be evident by all the one-line definitions. The conditional function (a variant of switch) is especially useful.
op.die <- function(op, d1, d2) {
if(missing(d2)) {
values <- op(d1$value)
probs <- d1$prob
} else {
values <- c(outer(d1$value, d2$value, FUN=op))
probs <- c(outer(d1$prob, d2$prob, FUN='*'))
}
die.normalize(values, probs)
}
"[.die" <- function(d1, i) sum(d1$prob[d1$value %in% i])
"==.die" <- function(d1, d2) op.die('==', d1, d2)
">.die" <- function(d1, d2) op.die('>', d1, d2)
"<=.die" <- function(d1, d2) op.die('<=', d1, d2)
"!.die" <- function(d) op.die(function(x) 1-x, d)
"+.die" <- function(d1, d2) op.die('+', d1, d2)
"-.die" <- function(d1, d2) op.die('-', d1, d2)
"*.die" <- function(d1, d2) op.die('*', d1, d2)
"/.die" <- function(d1, d2) op.die('/', d1, d2)
sign.die <- function(d) op.die(sign, d)
half <- function(d) op.die(function(x) floor(x/2), d)
conditional <- function(cond, dice) {
values <- unlist(sapply(dice, function(x) x$value))
probs <- unlist(sapply(1:length(cond$prob),
function(i) cond$prob[i] * dice[[i]]$prob))
die.normalize(values, probs)
}
(If one wanted to be efficient, which might be useful when working with large distributions, rep.die, +.die, and -.die could be specially rewritten to use convolutions. This is unlikely to be helpful in most applications, though, because the other operations would still need brute-force calculation.)
To enable study of the properties of distributions, here are some statistical summaries:
moment <- function(d, k) sum(d$value^k * d$prob)
mean.die <- function(d) moment(d, 1)
var.die <- function(d) moment(d, 2) - moment(d, 1)^2
sd.die <- function(d) sqrt(var.die(d))
min.die <- function(d) min(d$value)
max.die <- function(d) max(d$value)
As an example of their use, here is the health distribution for three damage rolls (the right hand plot in the first figure). The calculation of the total damage distribution is performed by x.3 <- health - rep(damage, 3) (pretty simple, right?) and the Normal approximation is computed via pnorm(x, mean.die(x.3), sd.die(x.3)).
plot(x.3 <- health - rep(damage, 3), type="b", xlim=l, lwd=2, xlab="Health",
main="After Three Hits")
curve(pnorm(x, mean.die(x.3), sd.die(x.3)), lwd=2, col="Red", add=TRUE)
abline(v=0, col="Gray")
abline(h = (x.3 <= nothing)[TRUE], col="Blue")
All this ought to port easily to C++.
Best Answer
If the outcome for each missile is independent of all other missiles (not necessarily realistic) and $p$ is the same for all missiles, then the number of missiles that hit is binomial$(N,p)$.
If $X$ is the number of hits then $$P(X=x)={N\choose x} p^x (1-p)^{N-x}\,,\quad x=0,1,2,...,N$$
See Wikipedia on the binomial distribution.
Similarly, under the same assumptions, the probability of at least one hit is 1-P(0 hits), and P(0 hits) is just $P(X=0)$ above.
To randomly generate the number of hits, you can use functions that randomly generate from that binomial distribution (I don't know what you have access to but libraries for this sort of thing are commonplace).
Alternatively, since partial hits are okay, you could randomly generate from a normal with mean $Np$, and variance $Np(1-p)$, and if the answer is below 0 call it 0, and if it's above N call it N. Nearly always you'll get a non-integer number in between.
You could then round to the nearest integer, or you can say "if it's within some distance $d$ of $x+\frac{1}{2}$, call it $x$ hits and one partial hit". So you might say d=0.25, for example, which means if you got a result of 5.2 you'd say 'five hits' but a result of 5.4 would be '5 hits and one partial hit', while 5.9 would be "six hits'. By playing with $d$ you can increase or decrease the frequency with which any partial hits show up. This will never give more than a single partial hit though.
If you don't want to assume independence, you'll have to specify the dependence you want.