What is the distribution of the peak time of the first hitting time process

Question

I need to find the distribution of the random variable $T_{peak}$ where $T_{peak}$ represents the peak time of the first hitting time process.

Detailed Explanation of the System: There are $N^{Tx}$ amount of emitted molecules from a specific point in 3D environment. Molecules diffuse in the environment according to the followings:

$$ r[t] = r[t-1] + (\Delta r_1, \Delta r_2, \Delta r_3)$$ $$ \Delta r_i \sim \mathcal{N}(0,\, 2D\Delta t)$$ where $r[t]$, $r_i$, $D$, and $\Delta t$ are the location vector at time $t$, $i$-th component of the location vector, diffusion coefficient, and the time step, respectively.

If there is an absorbing spherical trap at a distance $d$, the mean number of arriving/hitting molecules until time $t$ is:

$$ E[N^{absorb}(t)] = N^{Tx} \frac{r_{trap}}{d+r_{trap}} \, \text{erfc} \left( \frac{d}{\sqrt{4Dt}} \right) = N^{Tx} \frac{r_{trap}}{d+r_{trap}} \, 2\Phi \left( \frac{-d}{\sqrt{2Dt}} \right)$$ where $r_{trap}$ is the radius of the absorbing spherical trap.

When you plot $N^{absorb}(t)$ in small intervals, you get something like the following figure (Scaled Inverse Gaussian distribution)

And the expected value of the peak time of hitting time histogram is $$ E[T_{peak}] = \frac{d^2}{6D}$$

When simulating this diffusion process and the focusing on the absorbing times, $T_{peak}$ differs from simulation instance to instance, $T_{peak}$ is a random variable and I need to find the distribution of $T_{peak}$.

::Edit/New Text Start

If we consider time bins $\{t_0, t_1, ..., t_k, ...\}$ and denote the number of absorbed molecules in that bin as $N(t_k)=N^{absorb}(t_k^-)-N^{absorb}(t_k^+)$ where $t_k^-$ and $t_k^+$ are the left and the right end of the $k$-th time bin, we can define $T_{peak}$ as follows: $$T_{peak} = \arg\max\limits_{t_k} N(t_k)$$ with this definition, what is the distribution of $T_{peak}$. Is it Poisson or Binomial or Gaussian or Inverse Gaussian? Does it have well known structure and name?

Example Topology Figure:

::Edit/New Text End

P.S. These absorbed molecules are considered as the received signal in molecular communications and the peak time distribution of the received signal is important for many applications.

Answer 1

Description as multinomial distribution

You could view the distribution of the particles in the different time bins as a multinomial distribution where the mean amount of particles $\lambda_k$ in a particular time bin $t_k$ is $N(t_{k+1})-N(t_{k})$, or approximating this with the derivative of $N(t)$

$$\lambda_k = N^\prime(t) \Delta t = \frac{r_{trap}}{d+r_{trap}} \frac{n}{\sqrt{2Dt_k^3}} \phi\left(\frac{-d}{\sqrt{2Dt_k}}\right) \Delta t$$

where $n$ is the total number of particles.

---

Approximations

When the number of bins is large, when there is less correlation between pairs of bins, then you can approximate the multinomial distribution by individual binomial distributions.

When the number of particles is large then these binomial distributions can be approximated by poisson distributions or normal distributions.

---

Probability distribution for time bin with maximum number of molecules

With the distribution of number of molecules $n_k$ in the $k$-th time bin described by:

$$n_k \sim Pois(\lambda_k)$$

You can compute the probability for a particular bin to contain the extreme value as

$$P(\text{max at $t_k$}) = \sum_{i=0}^\infty \underbrace{f(i,\lambda_k)}_{\substack{\text{probability bin $k$} \\ \text{ contains $i$ molecules}}} \prod_{j \neq k} \underbrace{ F(i,\lambda_j) }_{\substack{\text{probability bin $j$} \\ \text{ contains at most $i$ molecules$$}}} $$

Where $f$ and $F$ are the density and cumulative distribution function of the Poisson distribution.

---

This is a bit of an horrible expression, but it can already be computed. I will see at a later time If I can make it look easier.

Computation example:

The image below shows the outcome of a simulation along with the theoretic distribution for the peak time (which seems to work well)

It resembles a curve that might be approximated with a normal curve (although it is not entirely the same).

# model parameters
dt <- 0.001
t <- seq(dt,0.3,dt)
n=700*1000
D = 1
d = 0.5

# model
ft <- n*d/sqrt(2*D*t^3)*dnorm(d/sqrt(2*D*t),0,1)
fmids <- n*d/sqrt(2*D*(t+dt/2)^3)*dnorm(d/sqrt(2*D*(t+dt/2)),0,1)
plot(t,ft*dt,type="l",lwd=1.5,lty=2)

# simulation
#    
# simulation by drawing from uniform distribution
# and converting to time by using quantile function of normal distribution
ps <- runif(n,0,1)                
ts <- 2*pnorm(-d/sqrt(2*D*t))     
sumn <- sapply(ts, FUN = function(tb) sum(ps < tb))
lines(t[-length(sumn)],sumn[-1]-sumn[-length(sumn)],col=4)

# compute t peak distribution
pmax = rep(0,length(t))         # probability of max
for (i in 1:3000) {
  Fj <- ppois(i,fmids*dt,log.p = TRUE)
  ss <- sum(Fj)
  for (j in 1:length(t)) {
    # doing the product by taking sum of log values
    pmax[j] <- pmax[j]+exp(dpois(i,fmids[j]*dt,log = TRUE) + ss - Fj[j])
  }
}

# compute whether it sums to 1
# you will see a bit higher because a maximum value may occur at multiple times
sum(pmax)

#multiple simulations
nsim <- 10000
histt <- c(0) #container with simulated maximums
pb <- txtProgressBar(title = "progress bar", min = 0,
                     max = nsim, style=3)
for (i in 1:nsim) {
  ps <- runif(n,0,1)
  ts <- 2*pnorm(-d/sqrt(2*D*t))
  sumn <- sapply(ts, FUN = function(tb) sum(ps < tb))
  histt <- c(histt,which(sumn[-1]-sumn[-length(sumn)] == 
                            max(sumn[-1]-sumn[-length(sumn)])))
  setTxtProgressBar(pb, i)
}
close(pb)
histt <- histt[-1]  # strip leading 0

h <- hist(histt*dt,breaks = t,xlim=c(0.025,0.06),main="peak time",
          xlab = "time")
lines(h$mids,pmax[-1]*nsim,col=2)
points(h$mids,pmax[-1]*nsim,col=2)

qqnorm(histt, pch=21,cex=0.5,col=1,bg=1)
qqline(histt)