Example 2: Bayesian model selection

In this example, designers of a multiserver system are unsure about the demand for service. Specifically, demand is either 'average' (typical scenario) with a number of arrivals pmf

{2: 0.3, 3: 0.7}

or 'high' with a number of arrivals pmf

{2: 0.2, 3: 0.3, 4: 0.5}

Based on experience, the prior probability of a 'high' demand scenario is initially set to 0.1. At time 5 the number in the system will be observed. If the posterior probability of the 'high' demand scenario is sufficiently high, management can begin the process of procuring more servers.

For purposes of this example, suppose that, at time 5, the observed value of the number of entities in the system is 12.

Calculating the posterior probabilities of the average and high demand scenarios proceeds in two steps:

Step 1. Run engine for five periods under each demand scenario and determine the respective probabilities of observing 12 entities in the system upon completion.
Step 2. Calculate posterior probabilities of each demand scenario using Bayes' theorem.

Here is the code to create the engine and specify its attributes:

from qplex import StandardMultiserver

engine = StandardMultiserver()
engine.service_duration_pmf = {1: 0.6, 2:0.3, 3: 0.1}
engine.number_of_servers = 2

avg_number_of_arrivals_pmf = {2: 0.3, 3: 0.7}
high_number_of_arrivals_pmf = {2: 0.2, 3: 0.3, 4: 0.5}
prior_prob_of_high = 0.1

Note that the number of arrivals pmf attribute is not initially assigned as there are two possibilities.

Here is the code to implement Step 1:

engine.mark()
engine.number_of_arrivals_pmf = avg_number_of_arrivals_pmf
for t in range(5):
    engine.step()
prob_outcome_avg = engine.get_number_of_entities_in_system_pmf()[12]

engine.restore()
engine.number_of_arrivals_pmf = high_number_of_arrivals_pmf
for t in range(5):
    engine.step()
prob_outcome_high = engine.get_number_of_entities_in_system_pmf()[12]

Here is the code to implement Step 2:

numerator = prob_outcome_high * prior_prob_of_high
denominator = prob_outcome_high * prior_prob_of_high + prob_outcome_avg * (1 - prior_prob_of_high)
posterior_prob_of_high = numerator/denominator
print('Posterior probability of high demand scenario = ', posterior_prob_of_high)

Here is the complete code:

from qplex import StandardMultiserver

engine = StandardMultiserver()
engine.service_duration_pmf = {1: 0.6, 2:0.3, 3: 0.1}
engine.number_of_servers = 2

avg_number_of_arrivals_pmf = {2: 0.3, 3: 0.7}
high_number_of_arrivals_pmf = {2: 0.2, 3: 0.3, 4: 0.5}
prior_prob_of_high = 0.1

engine.mark()
engine.number_of_arrivals_pmf = avg_number_of_arrivals_pmf
for t in range(5):
    engine.step()
prob_outcome_avg = engine.get_number_of_entities_in_system_pmf()[12]

engine.restore()
engine.number_of_arrivals_pmf = high_number_of_arrivals_pmf
for t in range(5):
    engine.step()
prob_outcome_high = engine.get_number_of_entities_in_system_pmf()[12]

numerator = prob_outcome_high * prior_prob_of_high
denominator = prob_outcome_high * prior_prob_of_high + prob_outcome_avg * (1 - prior_prob_of_high)
posterior_prob_of_high = numerator/denominator
print('Posterior probability of high demand scenario = ', posterior_prob_of_high)

And here is the output:

Posterior probability of high demand scenario = 0.6401739159505406

It is now more likely than not that the high demand scenario is the true scenario.