What is the interleaved probability like when two Gamma distribution processes fired together?

Question

As opposed to the similar question here: Is there a probabilistic (not analytical) argument for why the sum of independent Poissons is Poisson?

The difference is to consider the interval between two adjacent events follows the Gamma distribution given by two parameters $a$ and $b$, where the probability of having interval $D$ is $$ p(D|a;b)=\left(\frac{D^{a-1}}{\Gamma(a)}\right)\frac{e^{-D/b}}{b^a}, \ a,b>0, D\geq 0 $$ Note that there are now shape parameter $a$ and scale parameter $b$. As comparing to Poisson Distribution, we could get the estimated firing rate (number of events per second) is $\frac{1}{a*b}$ i.e. the intensity in the previous answer $\alpha = \frac{1}{a*b}$

Suppose similarly, the same signal is fired following two independent Gamma distributed processes $Gamma(a_1,b_1)$ and $Gamma(a_2,b_2)$ and the signal could not be distinguished from which process caused it. From the previous mentioned question's answer, when observing the firing rate of the signal, the mean rate would be equal to the sum of their rate $\frac{1}{a*b}=\frac{1}{a_1*b_1}+\frac{1}{a_2*b_2}$.

But, in this question, the Gamma Distribution has two separate parameters, could the parameters $a_1, b_1$ and $a_2, b_2$ be estimated separately based on the combined observation of measurements of intervals? Is the combined pattern still following a Gamma distribution?

In addition, intuitively it might doesn't matter whether the two processes' starting time aligns. Think about a special case, when $a_1 b_1=a_2 b_2$, but they starts with a tiny difference, comparing to another case 2 where they miss-align the starting time by $ab/2$, will this difference affect the resultant distribution of interval when $t \rightarrow \infty$ ? If the interleaved process keeps at case 1 as always in good alignment for the two internal sources generating the signal, the resultant distribution won't be a Gamma Distribution at all as there are two very probable peaks in the interval's PDF. It seems like a paradox that whether random process is random or following the statistical rule. But in nature, the probability of two internal sources have exact the same firing rate is zero. There will always be some drifts over time I guess ?

Answer 1

Question Description

Analysis

Use Case 1 as an example, named it as $D_{11}=t$

1. First, considering $D_{11}=t$ without the interference of Unit2, i.e., $D_{1}=t$ \begin{align} P(D_1=t) = PDF(t;\kappa_1,\theta_1) \end{align} 
2. Then, think about the latest Unit2 event. Assume that the latest event 2 occurs s second before the leading event, we get: \begin{align} P(D_2>t+s) &= \int_{t+s}^\infty PDF(\tau,\kappa_2,\theta_2) d\tau \\ &= 1 - CDF(t+s;\kappa_2,\theta_2) \end{align}
3. Since the two events are independent of each other, we can assume： \begin{align} s\sim U(0,\infty) \end{align} Finally, we have: \begin{align} P(D_{11}=t) &= Norm\bigg[ f(t) \bigg]\\ f_{11}(t)&=pdf_1(t)\int_0^{\infty}1-cdf_2(t+s)\ ds \end{align}

Similarly, we have: \begin{align} f_{11}(t)&=pdf_1(t)\int_0^{\infty}1-cdf_2(t+s)\ ds \\ f_{12}(t)&=(1-cdf_1(t))\int_0^{\infty}pdf_2(t+s)\ ds \\ f_{21}(t)&=(1-cdf_2(t))\int_0^{\infty}pdf_1(t+s)\ ds \\ f_{22}(t)&=pdf_2(t)\int_0^{\infty}1-cdf_1(t+s)\ ds \\ \end{align}

\begin{align} \\[3mm] P_{11}(t)&=\frac{pdf_1(t)\int_0^{\infty}1-cdf_2(t+s)\ ds }{a_1b_1+a_2b_2}\ \\[3mm] P_{12}(t)&=\frac{(1-cdf_1(t))\int_0^{\infty}pdf_2(t+s)\ ds}{a_1b_1+a_2b_2}\ \\[3mm] P_{21}(t)&=\frac{(1-cdf_2(t))\int_0^{\infty}pdf_1(t+s)\ ds}{a_1b_1+a_2b_2}\ \\[3mm] P_{22}(t)&=\frac{pdf_2(t)\int_0^{\infty}1-cdf_1(t+s)\ ds }{a_1b_1+a_2b_2}\ \\[3mm] \end{align}

Result

The statistical results on the left side are derived from five sets of Gamma distributions generated using MATLAB, showcasing the counts resulting of $[D_{11}, D_{12}; D_{21}, D_{22}]$from pairwise mixing scenarios. On the right side, the statistical outcomes are computed using the aforementioned formulas in Mathematica. It is evident that there is minimal disparity between the computed results and the statistical outcomes.

Discussion

This result is not perfectly. Ideally, it should be simplified into a more general form. Furthermore, it is desirable that the computed results do not require manual normalization and directly yield probabilities.

Answer 2

considering two conditions

To find the interleaved interval's distribution, we could consider four cases 11,12,21,22

We can also consider the two cases 1* and 2*.

• The waiting time untill the next point conditional on the current point being process I.
• The waiting time untill the next point conditional on the current point being process II.

Then mix the two together based on their relative frequency of occurrence.

Waiting time conditional on current event being a specific type.

If the current point is type I, then let's consider the waiting time distribution untill the next event I and the next evenent II. We call them respectively $f(t)$ and $g(t)$.

• The distribution $f(t)$ is the gamma distribution of the intervals in events of type I. $$f(t) = \text{pdf}_{gamma}(t,\alpha_1,\beta_2)$$

• The distribution $g(t)$ is as described in Distribution of conditional waiting time untill next event in a renewal process $$g(t) = \frac{\beta_2}{\alpha_2} (1-\text{cdf}_{gamma}(t,\alpha_2,\beta_2)) $$

The time untill the next event, either type I or type II, is the minimum of these two waiting times and is distributed as $$h(t) = f(t) (1-G(t)) + g(t) (1-F(t))$$

Worked out as code

The code below demonstrates the above formula in a concrete example.

### function to simulate the intervals 
sim = function(k1 = 1, r1 = 1, k2 = 1, r2 = 1, n = 10^6) {
  t1 =  rgamma(1,k1,r1)
  t2 =  rgamma(1,k2,r2)
  t = rep(NA,n)

  for (i in 1:n) {
     t[i] = min(t1,t2)
     if (t1<t2) {
       t1 = t1 + rgamma(1,k1,r1)
     } else {
        t2 = t2 + rgamma(1,k2,r2)
     } 
  }
  d = diff(t)
  return(d) 
}


### some settings
set.seed(1)

alpha1 = 6
beta1 = 30
lambda1 = beta1/alpha1    
alpha2 = 1
beta2 = 5
lambda2 = beta2/alpha2

### run simulate function 
x = sim(alpha1,beta1,alpha2,beta2)
### plot histogram 
hist(x,breaks = seq(0,max(x)+0.01,0.01), freq = 0, 
     main = "histogram of simulated interval times \n with superposed curve for approximation of distribution")

#######
## compute the distribution 
######
xs = seq(0,1,0.001)

f = dgamma(xs,alpha1,beta1)
F = pgamma(xs,alpha1,beta1)
g = beta2/alpha2*(1-pgamma(xs,alpha2,beta2))
G = 1+1/alpha2*(
   xs*beta2*(1-pgamma(xs,alpha2,beta2))-
     alpha2*(1-pgamma(xs,alpha2+1,beta2))
   )

h1 =    f*(1-G) + g*(1-F)

f = dgamma(xs,alpha2,beta2)
F = pgamma(xs,alpha2,beta2)
g = beta1/alpha1*(1-pgamma(xs,alpha1,beta1))
G = 1+1/alpha1*(
   xs*beta1*(1-pgamma(xs,alpha1,beta1))-
   alpha1*(1-pgamma(xs,alpha1+1,beta1))
   )

h2 =    f*(1-G) + g*(1-F)

h = (lambda1 * h1 + lambda2 * h2)/(lambda1+lambda2)

lines(xs,h, col = 2)