Example 3: Joint distributions

This example calculates the joint distribution of the number of entities in a system at times 1 and 3. For notational convenience, let Z1 and Z3 denote these random variables. Here is a high level outline of the logic:

• Step 1. Create the engine, advance it one step, mark to save the state, and obtain the distribution of Z1.

• Step 2. For each value of Z1, restore the internal state to its value at time 1 and apply the function apply_observation_value to condition the engine internal state on the value of Z1. Advance the engine two more steps to determine the conditional distribution of Z3 given the value of Z1. Then, for each value of Z3, print the joint probability P[Z1 = z1, Z3 = z3] = P[Z1 = z1] * P[Z3 = z3 | Z1 = z1].

Here is the code from Lesson 3 to create an engine and specify its attributes:

from qplex import StandardMultiserver

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

Here is the code from Lesson 4 to complete Step 1:

engine.step()
engine.mark()
pmf_1 = engine.get_number_of_entities_in_system_pmf()
print(pmf_1)

And here is the familiar output:

{2: 0.3000000000000001, 3: 0.7}

Here is the loop for Step 2:

for z1 in pmf_1:
    engine.restore()
    engine.apply_observation_value(z1)
    print(engine.get_number_of_entities_in_system_pmf())

When z1 is 2, the output is

{2: 1.0}

and when z1 is 3, the output is

{3: 1.0}

Here is the code to advance the engine to time 3:

engine.step()
engine.step()
pmf_3_given_1 = engine.get_number_of_entities_in_system_pmf()
print(pmf_3_given_1)

When z1 is 2, the output is

{2: 0.011663999999999999, 3: 0.08951236363636364, 4: 0.2600495795454546, 5: 0.35278537500000007, 6: 0.22281105681818172, 7: 0.05827762499999996, 8: 0.0048999999999999955}

and when z1 is 3, the output is shifted:

{3: 0.011663999999999999, 4: 0.08951236363636364, 5: 0.2600495795454546, 6: 0.35278537500000007, 7: 0.22281105681818172, 8: 0.05827762499999996, 9: 0.0048999999999999955}

Finally, here is the code to print the joint pmf:

for z3 in pmf_3_given_1:
    joint_prob = pmf_1[z1] * pmf_3_given_1[z3]
    print('P[ Z1 =', z1, ', Z3 =', z3, '] =', joint_prob)
print()

Here is the complete code:

from qplex import StandardMultiserver

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

engine.step()
engine.mark()
pmf_1 = engine.get_number_of_entities_in_system_pmf()

for z1 in pmf_1:
    engine.restore()
    engine.apply_observation_value(z1)

    engine.step()
    engine.step()
    pmf_3_given_1 = engine.get_number_of_entities_in_system_pmf()

    for z3 in pmf_3_given_1:
        joint_prob = pmf_1[z1] * pmf_3_given_1[z3]
        print('P[ Z1 =', z1, ', Z3 =', z3, '] =', joint_prob)
    print()

And here is the output:

P[ Z1 = 2 , Z3 = 2 ] = 0.003499200000000001
P[ Z1 = 2 , Z3 = 3 ] = 0.0268537090909091
P[ Z1 = 2 , Z3 = 4 ] = 0.07801487386363641
P[ Z1 = 2 , Z3 = 5 ] = 0.10583561250000005
P[ Z1 = 2 , Z3 = 6 ] = 0.06684331704545454
P[ Z1 = 2 , Z3 = 7 ] = 0.017483287499999993
P[ Z1 = 2 , Z3 = 8 ] = 0.001469999999999999

P[ Z1 = 3 , Z3 = 3 ] = 0.008164799999999998
P[ Z1 = 3 , Z3 = 4 ] = 0.06265865454545455
P[ Z1 = 3 , Z3 = 5 ] = 0.1820347056818182
P[ Z1 = 3 , Z3 = 6 ] = 0.24694976250000003
P[ Z1 = 3 , Z3 = 7 ] = 0.1559677397727272
P[ Z1 = 3 , Z3 = 8 ] = 0.040794337499999965
P[ Z1 = 3 , Z3 = 9 ] = 0.0034299999999999964