PDF for total duration of "on" state in a Poisson process with binary values

Question

I have a process where, at times given by a Poisson law of rate $\lambda$, a system picks two values - say "on" or "off" - with respective probabilities $p_0$ and $1 - p_0$. What is the PDF for the total duration of "on" state in a period $T \gg 1/\lambda$? I figured out that if the number of "on" states were fixed (say, to the average value $p_0 T \lambda$), it would be an Erlang distribution, since this would then be a sum of exponential laws for the durations - but this number actually fluctuates according to the Poisson distribution. Thank you for your help.,

Answer 1

If you start at time $0$ with the system picking an on/off value as it does later, the total on time $T_{on}$ has expectation $p_{on}T$.

There will always be a positive probability of the values $0$ and $T$ since there may be no new picks before $T$: specifically $\mathbb P(T_{on}=0)=(1-p_{on}) e^{-\lambda T}$ and $\mathbb P(T_{on}=T)=p_{on} e^{-\lambda T}$, with a continuous distribution between these values, so there is not going to be a simple solution for this mixed distribution.

Simulation suggests that its variance may be almost $2p_{on}(1-p_{on}) \frac{T}{\lambda}$ and the density and cdf curves may eventually not be far from normal distribution curves, remembering that $T_0$ is constrained to values in $[0,T]$ and has points of positive probability at the ends of this interval.

For example using R:

timeon <- function(T, lambda, pon){
  maxsamples <- ceiling(T*lambda +10*sqrt(T*lambda)+30)
  gaps <- rexp(maxsamples, rate=lambda)
  cumgaps <- cumsum(gaps)
  lastsample <- min(which(cumgaps > T))
  gaps <- gaps[1:lastsample]
  gaps[lastsample] <- gaps[lastsample]- (cumgaps[lastsample] - T)
  onoff <- sample(c(1,0), lastsample, replace=TRUE, prob=c(pon,1-pon))
  timeon <- sum(onoff * gaps) 
  return(timeon)
  }
    
T <- 2000
lambda <- 0.1
pon <- 0.3
set.seed(2023) # change if not trying to reproduce my numbers
sims <- replicate(10^5, timeon(T, lambda, pon))
mean(sims)
# 600.1637
T*pon # expectation
# 600
var(sims)
# 8375.701
2*pon*(1-pon)*T/lambda # possible function for variance
# 8400  

with a density (and normal approximation in red) looking like

plot(density(sims))
curve(dnorm(x, T*pon, sqrt(2*pon*(1-pon)*T/lambda)), add=TRUE, col="red")

and cdf (and normal approximation in red) looking like

plot.ecdf(sims)
curve(pnorm(x, T*pon, sqrt(2*pon*(1-pon)*T/lambda)), add=TRUE, col="red")