# Friendly probability puzzle solved with simulation

A friend posed to me a question to see if I could solve it. He says I solved it appropriately. I was hoping someone here could give me a review of the program I wrote to answer it. I figured it would be a good exercise in trying to do things as well as possible and try to follow conventions.

The problem

Aaron and Alexander are going to have a duel. Apparently, Alexander insulted Aaron and so Aaron has demanded a duel with pistols. The ground rules for the duel are that they will exchange single shots until someone is hit and thus the duel will end with the winner having been the first to hit the opponent. Both Aaron and Alexander have experience with duels and have tracked the frequency of their hits in past duels. freq Alexander hits = 0.6, freq Aaron hits = 0.5.

Using these frequencies as probabilities

A. What is the probability that Aaron will win the duel if he shoots first?

B. What is the probability that Aaron will win the duel if he shoots second?

After some digging and research, the problem can be solved treating the duel as an infinite geometric series and the result is A. 0.625 and B. 0.25. I'll save solving it that way in my code for another day. Instead, I solved it by simulating 10,000 duels using their hit probabilities. I know creating a class for a simple problem, but I figured it was a good exercise. Let me know what you think.

My code

import numpy as np

class Contestant:
"""Creates a contestant for a duel.

Attributes:
----------
name:   string
Name of the contestant

hit_prob:   NumPy array of float
Array containing the contestant's probability of a hit and probability of a miss

hit_arr:    NumPy array of int
Array containing the two possible outcomes of a shot, 1 = hit, 0 = miss.

shot_arr:   NumPy array of int
Array containing 100 predicted shots based on the hit_prob array

Methods:
-------
simulate_shots
Creates an array of predicted shot results for the shooter
"""

def __init__(self, name, hit_prob):
"""Initialize a contestant name, hit probability array, hit array, and an empty shot array."""

self.name = name
self.hit_prob = hit_prob
self.hit_arr = np.array([1, 0])
self.shot_arr = np.empty(100, dtype='int')

def simulate_shots(self, num_selection):
"""Class method to simulate a series of shots based on the contestant's hit probability.

Parameters:
----------
num_selection:  int
The number of shots to simulate.

Returns:
numpy.ndarray: Array of the simulated shots."""

# create an array of index numbers to select from
hit_indices = np.random.choice(len(self.hit_arr), size=num_selection, p=self.hit_prob)

# select those indices from the hit arrays
self.shot_arr = self.hit_arr[hit_indices]

def duel(first_shooter, second_shooter, volley):
"""Function to simulate a duel between the two contestants.

Parameters:
----------
first_shooter:  Contestant
Contestant that shoots first

second_shooter: Contestant
Contestant that shoots second

volley:   int
The current round of the duel.

Returns:
str: Name of the winner.
"""

if first_shooter.shot_arr[volley]:
return first_shooter.name

else:
if second_shooter.shot_arr[volley]:
return second_shooter.name

volley += 1

return duel(first_shooter, second_shooter, volley)

# the number of duels to simulate
number_of_duels = 10000

# probability arrays for the contestants
contestant_1 = Contestant('Aaron', np.array([0.5, 0.5]))
contestant_2 = Contestant('Alexander', np.array([0.6, 0.4]))

# the lists of wins
cont_1_shoots_first = []
cont_2_shoots_first = []

for duels in range(number_of_duels):

# simulate an array of shots for both contestants
contestant_1.simulate_shots(100)
contestant_2.simulate_shots(100)

# volley until there's a winner and append the winners name to the list
cont_1_shoots_first.append(duel(contestant_1, contestant_2, 0))
cont_2_shoots_first.append(duel(contestant_2, contestant_1, 0))

# calculate probabilities
print(f"If {contestant_1.name} shoots first, the probability he wins is "
f"{cont_1_shoots_first.count(contestant_1.name)/number_of_duels}.")

print(f"If {contestant_1.name} shoots second, the probability he wins is "
f"{cont_2_shoots_first.count(contestant_1.name)/number_of_duels}.")


Edit: "dual" to "duel"

Overall, it looks excellent:

• Good partitioning of code into class and functions
• Good layout and flow of the code
• Good naming of class, functions and variables
• Lint check (pylint) shows no major issues
• Great usage of docstrings/comments, especially on function inputs and return values

Your question text correctly uses the word "duel", and you even use it a couple times in the code comments. I think you really meant to use "duel" everywhere else in your code, instead of the word "dual". For example:

"""Creates a contestant for a dual.


You should mention the word "pistols" somewhere in your code comments, as you do in the question text.

Note: The question was changed after I posted this Answer. The original question used dual instead of duel in several places.

For the dual function, it would be helpful to mention in the comments that the function is called recursively. Also explicitly mention the ending condition for the recursion.

In the dual function, consider using an if/elif/else:

if first_shooter.shot_arr[volley]:
return first_shooter.name
elif second_shooter.shot_arr[volley]:
return second_shooter.name
else:
volley += 1
return duel(first_shooter, second_shooter, volley)


Firstly, it is a little more compact. Secondly, it makes it clear that one of three things will be returned.

I know creating a class for a simple problem, but I figured it was a good exercise.

I think what you are trying to say there is that you think it is overkill to use OOP here. Since you've already put in the work, it could allow for easier expansion of your code to add more features.

• Thank you so much for your comments. Yeah, I did mean the use of OOP. I originally had this all just in a single function. It wasn't anything that was going to be used anywhere. Just a fun problem a friend gave me but I thought it would be an opportunity to brush up and put a little work in. Yes, I did mean duel. in the dual function, is that particularly more readable for some? I did it the way I did to specifically make it more clear to me that first_shooter had to miss for the second_shooter.name to return.
– Dan
Feb 26 at 14:01
• @Dan: You're welcome. I updated my answer with a little more justification for the dual function. I find it a little easier to understand, but it is a matter of preference. Feb 26 at 14:09

ŃOT AN ANSWER, just a comment:

You do not need to solve 'infinite geometric series' once you enclose the geometric sequence into recursion.

Case A.

The first round, Aaron shoots.
– If he hits Alexander, he wins,
– but if he doesn't, Alexander shoots.
– – Then, if Alexander hits Aaron, Aaron loses,
– – but if Alexander misses, the next round starts and the situation is the same as at the start of the duel.

So Aaron can win either by winning in the first shot, or by missing the first shot and surviving the Alexander's response and winning the rest of duel.

Let $$\P_A\$$ – probability that Aaron wins,
$$\P_{A1}\$$ – probability that Aaron wins in the first round,
$$\Q_{A1}\$$ – probability that Alexander hits Aaron in his turn.

Then the probability that Aaron wins is

$$\P_A=P_{A1} + (1-P_{A1})\cdot (1-Q_{A1}) \cdot P_A \$$

so

$$\P_A (1 - (1-P_{A1})(1-Q_{A1})) = P_{A1}\$$

and, finally

$$\P_A = \frac {P_{A1}} { 1 - (1-P_{A1})(1-Q_{A1}) }\$$

where $$\ P_{A1} = 0.5, \ Q_{A1} = 0.6\$$.

Similarly for scenario B, you just need to multiply the above result by the probability of surviving the first Alexander's shot:

$$\ P_B = P_A\cdot (1-Q_{A1}) \$$

Often some simple maths saves us lots of encoding.

• Yep. I definitely understand the value in deriving the probability. I'm no mathematician, so I do not know what "enclosing in recursion" means, but I suspect that either I'm missing something or these are just different terms for the same thing. Or, maybe I'm mislabeling it, but if you consider that it could go 5 rounds, 6 rounds, 7 rounds, 8 rounds and that each no win round is represented by (1-Pa1)(1-Qa1), the equation you derived is that of an inf. geometric series. S = a/(1-r), where a = Pr(Aaron hit), and the common term r = Pr(Aaron miss) x Pr(Alex miss).
– Dan
Feb 26 at 15:07
• @Dan Precisely. And 'the recursion' hides in 'and the situation is the same as at the start of the duel', which manifests itself by $P_A$ re-appearing at the end of the first equation. Feb 27 at 6:49
• Is this a common thing? I saw another SE thread where someone was doing a similar dice problem and there was a similar answer (it may have been you that answered?). Basically, "You don't need an infinite geometric series", but when you looked at the equation derived, it was the same equation used for an infinite geometric series. Is it just that with probability, we call it something different?
– Dan
Feb 28 at 12:15

Your code is mostly well laid-out and clear. The only thing I'd suggest working on is naming, you use too many abbreviations and the duel method should have a verb in its name.

However, the logic needs some serious improvements.

# Don't define data twice

Probabilities must add to 1. Since there are only 2 outcomes for each shot, hit or miss, defining one of these probabilities implicitly defines the other. Numpy will raise an error if you try to work with invalid data, mitigating damage, but it's trivial to prevent invalid data in the first place.

# Don't simulate things that don't happen

You simulate 100 shots for each contestant (aside: this value of 100 is a magic number used multiple times. It should be a named constant instead). Running the simulation 2 million times shows that duels usually end after 1 to 9 rounds. This means that more than 90% of the data you generate is not used.

Also, if you want to initialize data before using a class, it's generally preferable to initialize it on instantiation or lazily when it's actually required, instead of requiring the caller to explicitly initialize it.

# Don't assume a random outcome will happen

It's extremely unlikely with your current input, but it's possible for the duel to last for more than 100 rounds, in which case your program will crash.

If you tweak the probabilities to lower the probability of a hit, this may become a real problem.

# Don't hang on to superfluous data

You save the outcome of every duel, while you only need to keep track of the total number of victories.

Repeatedly appending to lists is a relatively expensive operation, reallocating progressively more memory as the list grows. It's better to avoid it if you can.

# Putting it all together

Taking all this into account significantly improves code size and performance.

Smaller code is easier to maintain an extend, and easier to reason about to verify results.

Better performance means more simulations ran and better accuracy.

Here's my attempt:

from random import random

MAX_ITERATIONS = 1000

class Contestant:
'''
A constestant to a pistol duel
'''
def __init__(self, name, accuracy):
if accuracy < 0 or accuracy > 1:
raise ValueError('accuracy must be in range [0, 1]')
self.name = name
self.accuracy = accuracy

def simulate_duel(contestant_1, contestant_2, iteration=0):
'''
Simulate a pistol duel between two Constestant instances.
Returns the winner.
'''
if iteration > MAX_ITERATIONS:
raise OverflowError()
if random() < contestant_1.accuracy:
return contestant_1
return simulate_duel(contestant_2, contestant_1, iteration+1)

if __name__ == '__main__':
SIMULATION_COUNT = 1_000_000

contestant_1 = Contestant('Aaron', 0.5)
contestant_2 = Contestant('Alex', 0.6)

contestant_1_shoots_first = 0
contestant_2_shoots_first = 0

for _ in range(SIMULATION_COUNT):
if simulate_duel(contestant_1, contestant_2) is contestant_1:
contestant_1_shoots_first += 1
if simulate_duel(contestant_2, contestant_1) is contestant_1:
contestant_2_shoots_first += 1

print('Chance of Aaron winning if she goes first: ',
contestant_1_shoots_first / SIMULATION_COUNT)
print('Chance of Aaron winning if Alex goes first: ',
contestant_2_shoots_first / SIMULATION_COUNT)



Performance-wise, simulating 1 million duels with each contestant going first went down from over a minute to under a second.

• Thank you! I'll have to look into the calculation. I did that multiple times and found .625. Ill have to look at your code too and see where mine goes wrong because I still get a range around .625 with my simulation.
– Dan
Feb 28 at 11:08
• Ah! You have Aaron and Alex reversed. In my question and code Pr(Aaron hit) = 0.5 and Pr(Alex hit) = 0.6. The answer is 0.625 and 0.25...but your code is MUCH faster. I appreciate very much the input! The only reason I started out with defining both probabilities is I figured it would be faster to use NumPy to simulate a bunch of shots rather than do them one at a time, and using the NumPy library requires the entire prob space to be defined in that array. I wish I could mark both yours and @toolic responses as answers.
– Dan
Feb 28 at 12:22