¿Cuál es la distribución del tiempo pico del proceso de tiempo de primer impacto

Question

Necesito encontrar la distribución de la variable aleatoria $T_{peak}$ donde $T_{peak}$ representa el tiempo pico del proceso de primer tiempo de impacto.

Explicación detallada del sistema: Hay una cantidad $N^{Tx}$ de moléculas emitidas desde un punto específico en un entorno 3D. Las moléculas difunden en el entorno de acuerdo con lo siguiente:

$$ r[t] = r[t-1] + (\Delta r_1, \Delta r_2, \Delta r_3)$$ $$ \Delta r_i \sim \mathcal{N}(0,\, 2D\Delta t)$$ donde $r[t]$, $r_i$, $D$ y $\Delta t$ son el vector de ubicación en el tiempo $t$, el componente $i$-ésimo del vector de ubicación, el coeficiente de difusión y el paso de tiempo, respectivamente.

Si hay una trampa esférica absorbente a una distancia $d$, la cantidad media de moléculas que llegan/impactan hasta el tiempo $t$ es:

$$ 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)$$ donde $r_{trap}$ es el radio de la trampa esférica absorbente.

Cuando graficas $N^{absorb}(t)$ en intervalos cortos, obtienes algo así como la siguiente figura (Distribución inversa escalada gaussiana)

Y el valor esperado del tiempo pico del histograma de tiempo de impacto es $$ E[T_{peak}] = \frac{d^2}{6D}$$

Cuando se simula este proceso de difusión y se enfoca en los tiempos de absorción, $T_{peak}$ difiere de una instancia de simulación a otra, $T_{peak}$ es una variable aleatoria y necesito encontrar la distribución de $T_{peak}$.

::Edición/Nuevo texto Comienzo

Si consideramos los intervalos de tiempo $\{t_0, t_1, ..., t_k, ...\}$ y denotamos el número de moléculas absorbidas en ese intervalo como $N(t_k)=N^{absorb}(t_k^-)-N^{absorb}(t_k^+)$ donde $t_k^-$ y $t_k^+$ son el extremo izquierdo y el extremo derecho del intervalo de tiempo $k$-ésimo, podemos definir $T_{peak}$ de la siguiente manera: $$T_{peak} = \arg\max\limits_{t_k} N(t_k)$$ con esta definición, ¿cuál es la distribución de $T_{peak}$? ¿Es Poisson, Binomial, Gaussiana o Inversa Gaussiana? ¿Tiene una estructura y nombre bien conocidos?

Figura de topología de ejemplo:

::Edición/Nuevo texto Fin

P.D. Estas moléculas absorbidas se consideran como la señal recibida en comunicaciones moleculares y la distribución de tiempo pico de la señal recibida es importante para muchas aplicaciones.

Answer 1

Descripción como distribución multinomial

Podrías ver la distribución de las partículas en los diferentes intervalos de tiempo como una distribución multinomial donde la cantidad promedio de partículas $\lambda_k$ en un intervalo de tiempo particular $t_k$ es $N(t_{k+1})-N(t_{k})$, o aproximando esto con la derivada de $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$$

donde $n$ es el número total de partículas.

---

Aproximaciones

Cuando el número de intervalos es grande, cuando hay menos correlación entre pares de intervalos, entonces puedes aproximar la distribución multinomial por distribuciones binomiales individuales.

Cuando el número de partículas es grande, entonces estas distribuciones binomiales pueden aproximarse por distribuciones de Poisson o distribuciones normales.

---

Distribución de probabilidad para el intervalo de tiempo con el mayor número de moléculas

Con la distribución de número de moléculas $n_k$ en el $k$-ésimo intervalo de tiempo descrito por:

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

Puedes calcular la probabilidad de que un intervalo en particular contenga el valor extremo como

$$P(\text{máximo en $t_k$}) = \sum_{i=0}^\infty \underbrace{f(i,\lambda_k)}_{\substack{\text{probabilidad de que el intervalo $k$} \\ \text{contiene $i$ moléculas}}} \prod_{j \neq k} \underbrace{ F(i,\lambda_j) }_{\substack{\text{probabilidad de que el intervalo $j$} \\ \text{contiene a lo sumo $i$ moléculas$$}}} $$

Donde $f$ y $F$ son la función de densidad y la función de distribución acumulada de la distribución de Poisson.

---

Esta es una expresión un poco horrible, pero ya se puede calcular. Veré más adelante si puedo hacer que se vea más fácil.

Ejemplo de cálculo:

La imagen a continuación muestra el resultado de una simulación junto con la distribución teórica para el tiempo pico (que parece funcionar bien)

Se asemeja a una curva que podría aproximarse con una curva normal (aunque no es exactamente la misma).

# parámetros del modelo
dt <- 0.001
t <- seq(dt,0.3,dt)
n=700*1000
D = 1
d = 0.5

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

# simulación
#    
# simulación mediante generación de números aleatorios uniformes
# y convirtiéndolos a tiempo usando la función cuantil de la distribución normal
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)

# calcular distribución pico de t
pmax = rep(0,length(t))         # probabilidad del máximo
for (i in 1:3000) {
  Fj <- ppois(i,fmids*dt,log.p = TRUE)
  ss <- sum(Fj)
  for (j in 1:length(t)) {
    # haciendo el producto tomando la suma de los valores logarítmicos
    pmax[j] <- pmax[j]+exp(dpois(i,fmids[j]*dt,log = TRUE) + ss - Fj[j])
  }
}

# calcular si suma a 1
# se observará un poco más alto porque un valor máximo puede ocurrir en múltiples tiempos
sum(pmax)

# múltiples simulaciones
nsim <- 10000
histt <- c(0) #contenedor con máximos simulados
pb <- txtProgressBar(title = "barra de progreso", 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]  # eliminar el primer 0

h <- hist(histt*dt,breaks = t,xlim=c(0.025,0.06),main="tiempo pico",
          xlab = "tiempo")
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)