<jupyter_start><jupyter_text>Hyperexponential CaseThroughout this document, the following packages are required:<jupyter_code>import numpy as np
import scipy
import math
from scipy.stats import binom, erlang, poisson
from scipy.optimize import minimize
from functools import lru_cache<jupyter_output><empty_output><jupyter_text>Plot Phase-Type Fit<jupyter_code>from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
import matplotlib.pyplot as plt
def SCV_to_params(SCV):
    
    # weighted Erlang case
    if SCV <= 1:
        K = math.floor(1/SCV)
        p = ((K + 1) * SCV - math.sqrt((K + 1) * (1 - K * SCV))) / (SCV + 1)
        mu = K + (1 - p) * (K + 1)
    
        return K, p, mu
    
    # hyperexponential case
    else:
        p = 0.5 * (1 + np.sqrt((SCV - 1) / (SCV + 1)))
        mu = 1 # 1 / mean
        mu1 = 2 * p * mu
        mu2 = 2 * (1 - p) * mu
        
        return p, mu1, mu2
# for i in range(81):
#     SCV = 1 + 0.1 * i
#     print(round(SCV,2),SCV_to_params(SCV))
def density_WE(x, K, p, mu):
    return p * erlang.pdf(x, K, scale=1/mu) + (1 - p) * erlang.pdf(x, K+1, scale=1/mu)

def density_HE(x, p, mu1, mu2):
    return p * mu1 * np.exp(-mu1 * x) + (1 - p) * mu2 * np.exp(-mu2 * x)
x = np.linspace(0,4,1001)

def plot_f(SCV=1):
    
    if SCV <= 1:
        K, p, mu = SCV_to_params(SCV)
        f_x = density_WE(x, K, p, mu)
        title = f'SCV = {SCV}\n p = {p:.2f}, $K$ = {K}, $\mu$ = {mu:.2f}'
    else:
        p, mu1, mu2 = SCV_to_params(SCV)
        f_x = density_HE(x, p, mu1, mu2)
        title = f'SCV = {SCV}\n p = {p:.2f}, $\mu_1$ = {mu1:.2f}, $\mu_2$ = {mu2:.2f}'
    
    plt.plot(x,f_x)
    plt.title(title)
    plt.xlabel('$x$')
    plt.ylabel('density')
    plt.ylim(0,2)

interact(plot_f, SCV=(0.01,2,0.01));<jupyter_output><empty_output><jupyter_text>The recursion of the dynamic program is given as follows. For $i=1,\dots,n-1$, $k=1,\dots,i$, and $m\in\mathbb{N}_0$,\begin{align*}\xi_i(k,m) &= \inf_{t\in \mathbb{N}_0}\Big(\omega \bar{f}^{\circ}_{k,m\Delta}(t\Delta) + (1-\omega)\bar{h}^{\circ}_{k,m\Delta} +\sum_{\ell=2}^{k}\sum_{j=0}^{t}\bar{q}_{k\ell,mj}(t)\xi_{i+1}(\ell,j) +P^{\downarrow}_{k,m\Delta}(t\Delta)\xi_{i+1}(1,0) +P^{\uparrow}_{k,m\Delta}(t\Delta)\xi_{i+1}(k+1,m+t)\Big),\end{align*}whereas, for $k=1,\dots,n$ and $m\in \mathbb{N}_0$,\begin{align*}\xi_n(k,m) = (1-\omega)\bar{h}^{\circ}_{k,m\Delta}.\end{align*} We will implement this dynamic program step by step. First, we implement all functions in the equation above.Our formulas rely heavily on the survival function $\mathbb{P}(B>t)$ and $\gamma_z(t) = \mathbb{P}(Z_t = z\mid B>t)$:<jupyter_code>@lru_cache(maxsize=128)
def B_sf(t):
    """The survival function P(B > t)."""
    return p * np.exp(-mu1 * t) + (1 - p) * np.exp(-mu2 * t)

@lru_cache(maxsize=128)
def gamma(z, t):
    """Computes P(Z_t = z | B > t)."""
    
    gamma_circ = B_sf(t)
    
    if z == 1:
        return p * np.exp(-mu1 * t) / gamma_circ
    elif z == 2:
        return (1 - p) * np.exp(-mu2 * t) / gamma_circ<jupyter_output><empty_output><jupyter_text>Next, we implement $\bar{f}^{\circ}_{k,u}(t)$, which depends on $\bar{f}_{k,z}(t)$:<jupyter_code>@lru_cache(maxsize=128)
def f_bar(k,z,t):
    
    if z == 1:
        return sum([binom.pmf(m, k-1, p) * sigma(t, m+1, k-1-m) for m in range(k)])
    elif z == 2:
        return sum([binom.pmf(m, k-1, p) * sigma(t, m, k-m) for m in range(k)])

@lru_cache(maxsize=128)
def f_circ(k, u, t):
    return gamma(1, u) * f_bar(k, 1, t) + gamma(2, u) * f_bar(k, 2, t)<jupyter_output><empty_output><jupyter_text>In here, we need to evaluate the object $\sigma_{t}[m,k]$, which depends on $\rho_{t}[m,k]$:<jupyter_code>@lru_cache(maxsize=512)
def sigma(t,m,k):
    
    return (t - k / mu2) * erlang.cdf(t, m, scale=1/mu1) - (m / mu1) * erlang.cdf(t, m+1, mu1) + \
            (mu1 / mu2) * sum([(k-i) * rho_t(t, m-1, i) for i in range(k)])

@lru_cache(maxsize=512)
def rho_t(t,m,k):
    
    if not k:
        return np.exp(-mu2 * t) * (mu1 ** m) / ((mu1 - mu2) ** (m + 1)) * erlang.cdf(t, m+1, scale=1/(mu1 - mu2))
    elif not m:
        return np.exp(-mu1 * t) * (mu2 ** k) / ((mu1 - mu2) ** (k + 1)) * erlang.cdf(t, k+1, scale=1/(mu1 - mu2))
    else:
        return (mu1 * rho(t,a,m-1,k) - mu2 * rho(t,a,m,k-1)) / (mu1 - mu2)
    

@lru_cache(maxsize=512)
def rho(t,a,m,k):
    
    if not k:
        return np.exp(-mu2 * t) * (mu1 ** m) / ((mu1 - mu2) ** (m + 1)) * erlang.cdf(a, m+1, scale=1/(mu1 - mu2))
    elif not m:
        return np.exp(-mu1 * t) * (mu2 ** k) / ((mu1 - mu2) ** (k + 1)) * \
                    (erlang.cdf(t, k+1, scale=1/(mu1 - mu2)) - erlang.cdf(t-a, k+1, scale=1/(mu1 - mu2)))
    else:
        return (mu1 * rho(t,a,m-1,k) - mu2 * rho(t,a,m,k-1) - r(t,a,m,k)) / (mu1 - mu2)


@lru_cache(maxsize=512)
def r(t,s,m,k):
    return poisson.pmf(m,mu1*s) * poisson.pmf(k,t-s)<jupyter_output><empty_output><jupyter_text>We do the same for $\bar{h}^{\circ}_{k,u}(t)$, which only depends on $\bar{h}_{k,z}$:<jupyter_code>@lru_cache(maxsize=128)
def h_bar(k, z):

    if k == 1:
        return 0
    elif z <= K:
        return ((k - 1) * (K + 1 - p) + 1 - z) / mu
    elif z == K + 1:
        return ((k - 2) * (K + 1 - p) + 1) / mu

@lru_cache(maxsize=128)
def h_circ(k, u):
    return gamma(1, u) * h_bar() sum([gamma(z, u) * h_bar(k, z) for z in range(1, K+2)])<jupyter_output><empty_output><jupyter_text>The next objective is to implement $\bar{q}_{k\ell,mj}(t)$. This function depends on $q_{k\ell,z,v}(t)$, which depends on $\psi_{vt}[k,\ell]$: TODO<jupyter_code># TODO
poisson.pmf(3,0)<jupyter_output><empty_output><jupyter_text>Finally, we implement the remaining transition probabilities $P^{\uparrow}_{k,u}(t)$ and $P^{\downarrow}_{k,u}(t)$:<jupyter_code># @lru_cache(maxsize=128)
def P_up(k, u, t):
    """Computes P(N_t- = k | N_0 = k, B_0 = u)."""
    return B_sf(u + t) / B_sf(u)

@lru_cache(maxsize=128)
def P_down(k, u, t):
    """Computes P(N_t- = 0 | N_0 = k, B_0 = u)."""
    return sum([binom.pmf(m, k, p) * Psi(t, m, k-m) for m in range(k+1)])

@lru_cache(maxsize=128)
def Psi(t, m, k):
    return erlang.cdf(t, m, scale=1/mu1) - mu1 * sum([rho_t(t, m-1, i) for i in range(k)])
erlang.cdf(0,1,1)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.21321417172039578, 'fraction_numerical': 0.03702521544845196, 'mean_word_length': 3.0694805194805195, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 37, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '668208', 'license': 'MIT', 'n_tokens_mistral': 2944, 'n_tokens_neox': 2836, 'n_words': 624, 'path': '.ipynb_checkpoints/Hyperexponential Case-checkpoint.ipynb', 'repo_name': 'Roshanmahes/Appointment-Scheduling', 'size': '6266'}