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

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 ?

• Note that this question is also different to stats.stackexchange.com/questions/192067/… Because that question is assuming samples were drawn from two processes, their values only depend on how much portion of one process dominates. But here, the whole thing gets mixed and causing the observed sample value (i.e. the intervals) to be different. It is not putting all samples together in a bag after they were drawn. Each event might cause the interval observed to be altered or not depends on the two events source Commented Apr 19 at 2:05
• To find the interleaved interval's distribution, we could consider four cases 11,12,21,22 indicating the leading and ending edges source. P(1x) and P(2x) should be E1/(E1+E2) where E1=1/(a1*b1) i.e. intensity. And in 11 and 22 cases, the interval should follow original Gamma(a1,b1) and Gamma(a2,b2). But I have no idea about the interval in cases 12 and 21. Commented Apr 19 at 3:58
• I would think that the interarrival times of the combined process are neither gamma distributed nor independent. But you can probably do Bayesian inference of all four parameters by introducing some latent variables indicating whether each arrival time belongs to the first or second underlying process. You can then update these indicator variables along with the four main model parameters using Gibbs sampling to obtain samples from the joint posterior. Commented Apr 20 at 18:23
• Indeed, the superposition of two renewal processes (which is what you've got here) is not another renewal process, see en.wikipedia.org/wiki/Renewal_theory Commented Apr 20 at 20:02
• @JarleTufto - you might want to consider expanding your two comments into an answer, since they collectively do answer the questions! Commented Apr 20 at 21:09

# 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.

• In the derivation of $P(D_2>t)$ you consider two parts: (1) the conditional situation that the last event occured $s$ seconds ago $P(D_2>t+s)$ and (2) the probability density of the last event occuring $s$ seconds ago. In that derivation you suppose (1) that $P(D_2>t+s)$ is like the original interval time (but the interval is longer instead) and (2) that the distribution of $s$ is an improper uniform distribution. The end result is correct, but those two elements together are not correct; the errors happen to cancel each other out, but it is not a right motivation. Commented Apr 24 at 11:53

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