1
$\begingroup$

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()

enter image description here

$\endgroup$
0
$\begingroup$

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()

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.