The basic rules of engagement provide the probability distribution for transitions from $a$ attacking armies and $d$ defending armies to $a^\prime$ attackers and $d^\prime$ defenders, where $0 \le a^\prime \le a$ and $0 \le d^\prime \le d$. Beginning with $A$ attacking armies and $D$ defending armies, there are therefore $(A+1)(D+1)$ possible states indexed by $a=0, 1, \ldots, A$ and $d=0, 1, \ldots, D$. This is a Markov Chain on these states. All the information about outcomes is therefore contained in its transition matrix $\mathbb{P}_{A,D}$. In particular, after $k$ transitions ("turns" in the battle) from an initial state vector $s$, the distribution of resulting states $s^\prime$ is given by $\mathbb{P}_{A,D}^k\cdot s$. Because $\mathbb{P}_{A,D}$ can be large but $k$ is likely small, it may be most expedient simply to calculate this distribution iteratively via
$$s(k) = \mathbb{P}_{A,D}\cdot s(k-1)$$
beginning with the initial state $s(0) = s$.
These ideas were used to produce graphs of the winning chances for the attacker (blue) and defender (red) as a function of the number of turns played. (The gray lines are the chances that the battle remains unresolved after each turn.) The total time needed to do these computations using the brute-force method described below was five seconds. This is more than competitive with a Monte-Carlo simulation. For much larger numbers of armies, simulation will be superior in its use of time and RAM, but even so it is helpful to have a way to compute exact results, if only to test a simulator.
It is interesting that the attacker needs relatively few additional armies in order to secure a huge advantage. In particular, one might naively expect to need $0.7/0.6$ times as many attacking armies. The randomness of the outcomes magnifies a small initial superiority, however. Also, attackers have one slight advantage: they still win if they are wiped out, provided the defenders are also wiped out. (The actual rules might differ on this account; if so, the code would need slight modification to accommodate any additional restrictions concerning minimum numbers of attackers.)
An Algorithm
A solution with brute-force computation takes only a few seconds even when $A=100$, $D=60$, and $k=5$ as in the question. The only difficulties concern how to index the states. Because R
appears to be preferred and it organizes arrays by columns, I have written the transition vectors as columns and multiply state vectors by the transition matrix from the left (which is the opposite of the usual convention). I index the states $1, 2, \ldots, (A+1)(D+1)$ by letting the first coordinate vary fastest from largest to smallest. For instance, with $(A,D)=(2,1)$ the states are $(2,1), (1,1), (0,1), (2,0), (1,0), (0,0)$. The transition matrix $\mathbb{P}_{2,1}$ is computed in R
as
From
To 2:1 1:1 0:1 2:0 1:0 0:0
2:1 0.048 0.00 0 0 0 0
1:1 0.112 0.12 0 0 0 0
0:1 0.000 0.28 1 0 0 0
2:0 0.252 0.00 0 1 0 0
1:0 0.588 0.18 0 0 1 0
0:0 0.000 0.42 0 0 0 1
using the command mc(2, 1)$transition
(whose source is given below). For instance, the entry in row 1:0
and column 2:1
asserts there is a chance of 0.588
that two armies attacking one defender (the state $(2,1)$) will, in one turn, reduce to one attacking army and no defending armies (the state $(1,0)$).
To check the validity of this calculation, we need to know the rules. The transition from $(a,d)$ yields a pair of random variables $(a^\prime, d^\prime)$. Let $X$ be a Binomial$(a, p)$ variable, where by default $p = 0.6$, and $Y$ be an independent Binomial$(d, q)$ variable with $q = 0.7$ by default. Then $a^\prime = \max(0,a-Y)$ and $d^\prime = \max(0,d-X)$. We say that the attacker "wins" provided $a\gt 0$ and $d^\prime=0$ and that the attacker "loses" provided $a\gt 0,$ $a^\prime=0$, and $d^\prime\gt 0$. In the case of a win or loss the engagement ends. Otherwise it may continue (at the attacker's option).
Returning to the example, beginning in the state $(2,1)$ we find that $X$ has a Binomial$(2,0.6)$ distribution, whence
$$\Pr(X=0)=0.4^2=0.16,\Pr(X=1)=2(0.4)(0.6)=0.48, \Pr(X=2)=0.6^2=0.36.$$
Consequently
When $X=0, d^\prime=\max(0,1-0) = 1$, so ${\Pr}_{2,1}(d^\prime=1) = 0.16.$
When $X=1, d^\prime=\max(0,0)=0$; and
When $X=2, d^\prime=\max(0,-1)=0$, so ${\Pr}_{2,1}(d^\prime=0) = 0.48+0.36=0.84.$
Likewise $Y$ has a Binomial$(1,0.7)$ distribution, giving it a $0.3$ chance of being $0$ and a $0.7$ chance of being $1$. Accordingly
${\Pr}_{2,1}(a^\prime=2)=\Pr(Y=0)=0.3$ and
${\Pr}_{2,1}(a^\prime=1)=\Pr(Y=1)=0.7$.
Because $X$ and $Y$ are independent, the chances multiply. For example, the chance of arriving at $(a^\prime, d^\prime) = (1,0)$ is
$${\Pr}_{2,1}\left((a^\prime, d^\prime)\right) = {\Pr}_{2,1}(a^\prime=1){\Pr}_{2,1}(d^\prime=0) = 0.7\times 0.84 = 0.588.$$
The other entries in the transition matrix are similarly verified. Notice that the terminal states $(0,1), (2,0), (1,0),$ and $(0,0)$ just make transitions to themselves. The first is a loss for the attacker while the rest are counted as wins.
Working code
To implement this I have used auxiliary functions fight
to compute the probabilities, transition
to produce the distribution of transitions from a single state, and mc
to assemble those into a transition matrix. The function battle
repeatedly applies this transition matrix to the initial distribution where state $(A,D)$ has probability $1$ and all other states have probability $0$. It summarizes the attacker's chances of winning and losing. Finally, by iterating battle
over plausible choices of $A$, it becomes a tool to help with game strategy, which amounts to selecting an appropriate number of armies from among those available in order to attack a defender. The only applicable rule is that the number of armies must be strictly less than the number of available (because any attacking armies available after the engagement will be moved into the defender's territory but at least one army must be left behind).
(This modular implementation makes it relatively easy to check the calculations by hand, at least for small cases. Doing so increases the confidence that the calculations are correct. Moreover, by limiting each module (R
function) to less than a half dozen executable lines, each is sufficiently simple to undergo careful testing and review. I admit my testing has been cursory, but it did include a wide range of values of $A$, $D$, $p$, and $q$ as well as spot-checks like the example given above.)
The output of this code consists of an array of graphs showing the winning and losing chances of the attacker as a function of the number of turns. It also reports the total amount of time needed to do the calculations for each graph. On this workstation it takes under 2.5 seconds to perform the calculations for $A=100$ attackers, $D=60$ defenders, and $k=5$ turns (using the command battle(100, 60, 5)
). Easy modifications of battle
will enable deeper analysis, such as evaluating the numbers of attacking armies that are likely to survive any engagement: such information, more than the mere chance of winning after $k$ turns, is more important to developing an optimal game strategy.
fight <- function(n.attack, n.defend, n.defend.max, attack.p) {
#
# Returns a probability vector for the number of surviving defenders
# indexed by n.defend.max down to 0.
# NB: This is the only part that needs modification if a "luck factor"
# is introduced.
#
p.casualty <- dbinom(0:n.attack, n.attack, attack.p)
n.survive <- n.defend - (0:n.attack)
if (n.attack > n.defend) {
# Collapse the negative values into a single state with no survivors.
p.casualty <- c(p.casualty[n.survive > 0], sum(p.casualty[n.survive <= 0]))
n.survive <- n.survive[n.survive >= 0]
}
# Pad the return vector with zeros, fore and aft, as needed.
c(rep(0, n.defend.max-n.defend), p.casualty, rep(0, max(0, n.defend-n.attack)))
}
transition <- function(attack, defend, attack.max, defend.max, attack.p, defend.p) {
#
# Returns the transition probabilities from the state (attack, defend),
# in descending order with `attack` changing fastest.
#
a <- fight(defend, attack, attack.max, defend.p)
d <- fight(attack, defend, defend.max, attack.p)
return(as.vector(outer(a, d)))
}
mc <- function(attack.max, defend.max, attack.p=0.6, defend.p=0.7) {
#
# Returns the transition matrix for a round in which attack.max armies
# (not including any reserved ones) engage defend.max armies.
# Transitions are in *columns* (not rows, which is the convention).
# The matrix is square with dimensions (attack.max+1)(defend.max+1).
#
# Also returns indicator vectors $wins and $loses showing which
# states correspond to wins and losses for the attacker, respectively.
#
i <- expand.grid(A=attack.max:0, D=defend.max:0) # All states
x <- apply(i, 1, function(ad, ...) transition(ad[1], ad[2], ...),
attack.max=attack.max, defend.max=defend.max,
attack.p=attack.p, defend.p=defend.p)
# Name the indexes in `x` to assist human reading of the results.
s <- paste(i$A, i$D, sep=":")
dimnames(x) <- list(To=s, From=s)
return(list(transition=x, n=dim(i)[1], wins=(i$D==0), loses=(i$A==0 & i$D > 0)))
}
battle <- function(attack, defend, k=5, ...) {
#
# Conduct a battle of `attack` attacking armies against `defend` defending
# armies for `k` turns.
# Return an array whose rows (indexed by `k`) give chances of the attacker
# winning and losing at that turn.
#
# (The code is readily modified to report the distributions of numbers of
# remaining armies.)
#
if (attack <= 0 || defend <= 0) stop("Both army counts must be positive.")
#
# Find the probability distributions after 1, 2, ..., k turns.
#
x <- mc(attack, defend, ...) # Transition matrix structure
y <- c(1, rep(0, x$n-1)) # $Beginning distribution: all is in state 1.
p <- matrix(NA, k, 3,
dimnames=list(Turns=1:k, Outcome=c("Wins", "Loses", "Undecided")))
for (i in 1:k) {
y <- x$transition %*% y # $
# Summarize this turn.
p[i, "Wins"] <- sum(y[x$wins]) # $
p[i, "Loses"] <- sum(y[x$loses]) # $
p[i, "Undecided"] <- 1 - (p[i, "Wins"] + p[i, "Loses"])
}
return(p)
}
#
# Study near-equal battles involving a given number of defenders.
#
d <- 60
k <- 5
p.a <- 0.6; p.d <- 0.7
times <- numeric(0)
par(mfrow=c(2,2))
for (a in round(d * exp(seq(log(63/60), log(75/60), length.out=4)))) {
u <- system.time(p <- zapsmall(battle(a, d, k, attack.p=p.a, defend.p=p.d)))
times <- c(times, u[3])
plot(p[, "Wins"], type="l", ylim=c(0,1), col="Blue",
xlab="Turns", ylab="Probability",
main=paste(a, "Attackers vs.", d, "Defenders"),
sub=paste("p(Attack) =", p.a, "p(Defense) =", p.d))
lines(p[, "Undecided"], col="Gray")
lines(p[, "Loses"], col="Red")
}
times
p=0.6
when attacking andp=0.7
when defending. I now understand that from the way I give the example it might seem like it depends on troops, but it was just posed like this for illustrative purposes. $\endgroup$