# Unsure whether my continuous time Markov Chain distribution is correct

I' am reading Introduction to Stochastic Processes by Lawler and have hit problem 3.3(c) that I' am not sure I have correct.

$\textbf{3.3}$ Suppose $X_t$ and $Y_t$ are independent Poisson processes with parameters $\lambda_1$ and $\lambda_2$, respectively, measuring the number of calls arriving at two different phones. Let $Z_t=X_t+Y_t$.

$\textbf{(c)}$ Let T denote the first time that at least one call has come from each of the two phones. Find the density and distribution function of the random variable T.

So here is my answer for the density function. Let $X_1$ and $Y_1$ denote the first time a call as been received for each process. Then the density function is, \begin{align} P(T=t)&=P(T_{X_1}=t,T_{Y_1}\leq t) + P(T_{Y_1}=t,T_{X_1}\leq t) \\ &=P(T_{X_1}=t)(1-P(T_{Y_1}>t)+P(T_{Y_1}=t)(1-P(T_{X_1}>t) \\ &=\lambda_1e^{-\lambda_1t}(1-e^{-\lambda_2t})+\lambda_2e^{-\lambda_2t}(1-e^{-\lambda_1t}) \\ &=\lambda_1e^{-\lambda_1t}+\lambda_2e^{-\lambda_2t}-e^{-(\lambda_1+\lambda_2)t}(\lambda_1 + \lambda_2) \end{align}

I've done some simulations, but it does not seem to fully support the PDF. Here is the code and image:

from collections import Counter
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats

N=50000
legend=[]
for lambda1 in np.arange(0.5,2,0.5):
for lambda2 in np.arange(0.5,2,0.5):
legend.append('lambda1='+str(lambda1)+',lambda2='+str(lambda2))
lambda1 = float(lambda1)
lambda2 = float(lambda2)

X1 = stats.expon.rvs(scale=1.0/lambda1,size=N)
X2 = stats.expon.rvs(scale=1.0/lambda2,size=N)
X = np.column_stack((X1,X2)).T

density = []

for t in np.arange(0.5,10,0.1):
count = ((X>=(t-0.5)) & (X<t)).sum(axis=0)
count_1 = X[:,count>=1]
count_2 = sum((count_1<t).sum(axis=0)==2)
density.append(count_2/float(N))

density_predicted = []
for t in np.arange(0.5,10,0.1):
temp = lambda1 * np.e**(-lambda1*t) + lambda2*np.e**(-lambda2*t) - np.e**(-(lambda1 + lambda2)*t)*(lambda1 + lambda2)
density_predicted.append(temp)

plt.plot(density,density_predicted)
plt.xlim(0,1)
plt.ylim(0,1)
plt.legend(legend,prop={'size':7})
plt.xlabel('Simulated density')
plt.ylabel('Predicted density')
plt.savefig('3.3.c.png')
plt.close()
plt.clf()


OK I think I solved it. My mistake was in simulating the pdf. I was not computing a density, but rather a proportion of events that occurred within a time window $\triangle t$. I just needed to divide that proportion by $\triangle t$. See the amended code below:

import matplotlib.pyplot as plt
import numpy as np
from scipy import stats

lambda_interval = 0.5
N=50000
legend=[]
delta_t = 0.1
for lambda1 in np.arange(2,2.5,lambda_interval):
for lambda2 in np.arange(2,2.5,lambda_interval):
legend.append('lambda1='+str(lambda1)+',lambda2='+str(lambda2))
lambda1 = float(lambda1)
lambda2 = float(lambda2)

X1 = stats.expon.rvs(scale=1.0/lambda1,size=N)
X2 = stats.expon.rvs(scale=1.0/lambda2,size=N)
X = np.column_stack((X1,X2)).T

density = []

for t in np.arange(0.5,10,delta_t):
count = ((X>=(t-delta_t)) & (X<t)).sum(axis=0)
count_1 = X[:,count>=1]
count_2 = sum((count_1<t).sum(axis=0)==2)
density.append(count_2/float(N)/delta_t)

density_predicted = []
for t in np.arange(0.5,10,delta_t):
temp = lambda1 * np.e**(-lambda1*t) + lambda2*np.e**(-lambda2*t) - np.e**(-(lambda1 + lambda2)*t)*(lambda1 + lambda2)
density_predicted.append(temp)

plt.plot(density,density_predicted)
plt.legend(legend,prop={'size':7})
plt.xlabel('Simulated density')
plt.ylabel('Predicted density')
plt.savefig('3.3.c.png')
plt.close()
plt.clf()