text_plain_output_1.png image_output_1.png
import matplotlib.pyplot as plt

import numpy as np

import numpy as np
import numpy.linalg
from sklearn import linear_model
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
import matplotlib.pyplot as plt
def theta_m2_m1(n2=1, n1=1, a=1, b=1, c=1, mu_1=1, mu_2=1):
    p = (a / mu_1) ** n1 * ((a + b + c) / (mu_1 + mu_2)) ** n2
    return p
def theta_m1_m2(n2=1, n1=1, a=1, b=1, c=1, mu_1=1, mu_2=1):
    p = (b / mu_2) ** n2 * ((a + b + c) / (mu_1 + mu_2)) ** n1
    return p
def theta_m2_0m1(n2=1, a=1, b=1, c=1, mu_1=1, mu_2=1):
    p = mu_1 / mu_2 * ((a + b + c) / (mu_1 + mu_2)) ** n2
    return p
def theta_m1_0m2(n1=1, a=1, b=1, c=1, mu_1=1, mu_2=1):
    p = mu_1 / mu_2 * ((a + b + c) / (mu_1 + mu_2)) ** n1
    return p
def theta_m2(n2=1, a=1, b=1, c=1, mu_1=1, mu_2=1):
    p = (mu_1 + mu_2) / (a + c) * (b / mu_2) ** n2
    return p
def theta_m1(n1=1, a=1, b=1, c=1, mu_1=1, mu_2=1):
    p = (mu_1 + mu_2) / (b + c) * (a / mu_1) ** n1
    return p
def theta_0s(a=1, b=1, c=1, mu_1=1, mu_2=1):
    A = np.array([[0, a + b + c + mu_2, -(b + c / 2)], [a + b + c + mu_1, 0, -(a + c / 2)], [-mu_1, -mu_2, a + b + c]])
    b = np.array([[(mu_1 + mu_2) * b / (a + c) + 2 * mu_1 ** 2 / mu_2], [(mu_1 + mu_2) * a / (b + c) + 2 * mu_1], [0]])
    A_inv = np.linalg.inv(A)
    ps = np.matmul(A_inv, b)
    return ps
def calc_visschers(a=1, b=1, c=1, mu_1=1, mu_2=1, k=100):
    p_list = []
    t0 = theta_0s(a, b, c, mu_1, mu_2)[2][0]
    p_list.append(t0)
    t1busy = theta_0s(a, b, c, mu_1, mu_2)[0][0]
    p_list.append(t1busy)
    t2busy = theta_0s(a, b, c, mu_1, mu_2)[1][0]
    p_list.append(t2busy)
    tbothbusy = theta_m2_0m1(n2=0, a=a, b=b, c=c, mu_1=mu_1, mu_2=mu_2) + theta_m1_0m2(n1=0, a=a, b=b, c=c, mu_1=mu_1, mu_2=mu_2)
    p_list.append(tbothbusy)
    for n in range(1, k):
        for m in range(1, k):
            t1 = theta_m2_m1(n2=m, n1=n, a=a, b=b, c=c, mu_1=mu_1, mu_2=mu_2)
            t2 = theta_m1_m2(n2=m, n1=n, a=a, b=b, c=c, mu_1=mu_1, mu_2=mu_2)
            p_list.append(t1)
            p_list.append(t2)
        t3 = theta_m2_0m1(n2=n, a=a, b=b, c=c, mu_1=mu_1, mu_2=mu_2)
        p_list.append(t3)
        t4 = theta_m1_0m2(n1=n, a=a, b=b, c=c, mu_1=mu_1, mu_2=mu_2)
        p_list.append(t4)
        t5 = theta_m2(n2=n, a=a, b=b, c=c, mu_1=mu_1, mu_2=mu_2)
        p_list.append(t5)
        t6 = theta_m1(n1=n, a=a, b=b, c=c, mu_1=mu_1, mu_2=mu_2)
        p_list.append(t6)
    x = np.array(p_list)
    x = x / x.sum()
    probs_array = x[0:4]
    return x

rlist = []
for k in range(100, 3000, 10):
    r = calc_visschers(a=4.75, b=4.75, c=0, mu_1=10, mu_2=10, k=k)
    rlist.append(r[0])
plt.plot(range(100, 3000, 10), rlist)