10
$\begingroup$

Suppose that $X$ measures the half-life of a radioactive element, with decay rate $\lambda$ (per unit of time).

Starting from a population of $N$ particle, I believe you can model the number of particles that decay over a unit of time as a $Poisson(N\lambda)$, and that the time until the first decay is detected is modeled using an $Exponential(N\lambda)$.

My question is: How can we model the distribution of the half-life $X$, i.e., starting from a population of $N$ particles, the amount of time it takes for $N/2$ of the particles to have decayed ?

Assuming the decays are independent from one another, it seems an Erlang distribution could be appropriate https://en.wikipedia.org/wiki/Erlang_distribution? Should I then use $f(x,N/2,N\lambda)$ (where the notation refers to the Wikipedia article)? If the half-life intuition is correct, it would seem $f(x,N/2,N\lambda)$ should then be invariant to the choice of $N$, but that's not obvious to me looking at the algebraic for of the Erlang's pdf. Is that correct? If so, there should be a single distribution $g = f(x,N/2,N\lambda)$ for all $N$ that models half-lifes?

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ I am not sure you are correct that the number of particles that decay over a unit of time has a Poisson distribution, partly because a Poisson distribution has unbounded support but you cannot have more than all the particles decay. $\endgroup$
    – Henry
    Commented Apr 22, 2022 at 15:29
  • $\begingroup$ @Henry: Good point. In practice, for this kind of application, $N$ is typically very large (number of particles in a batch of matter) and $\lambda$ relatively small. So I am hoping it would still provide a good approximation? Independence issues aside, to be more precise, I guess I could use a $Bi(N, \lambda)$? $\endgroup$
    – FZS
    Commented Apr 22, 2022 at 16:08
  • $\begingroup$ If $N$ is large and $\lambda$ is small then it will be a long time until half of the particles have decayed. Initially you will be seeing about $N \lambda$ decays per unit time and close to a Poisson process. The issue is that until this rate of decays starts to slow down noticeably, you will have no good way of estimating the overall half-life of about $\frac{\log_e(2)}{\lambda}$, and once it does slow down it will no longer be a homogeneous Poisson process. $\endgroup$
    – Henry
    Commented Apr 22, 2022 at 16:35
  • 1
    $\begingroup$ Hmm... I wonder, then, how physicists have managed to measure the half-lives of long-lived isotopes like $U^{238}$ (half life, $4.5\times 10^9$ years) with such accuracy. Answer: any macroscopic mass contains a huge number of particles so even over a short period of time a very large number of decays occur, providing precise estimation of $\lambda.$ @Henry $\endgroup$
    – whuber ♦
    Commented Apr 22, 2022 at 16:39
  • $\begingroup$ @whuber $N\lambda$ is easy to measure reasonably accurately if you have enough decays: just count the number of decays in a second or a day or whatever. The question is separating $N$ from $\lambda$ since you are interested in $\lambda$. For long-lived isotopes, you create a pure sample with a known total mass and a known atomic mass, dividing one by the other and multiplying by Avogadro's number then gives you $N$. $\endgroup$
    – Henry
    Commented Apr 22, 2022 at 16:55

1 Answer 1

Reset to default
14
$\begingroup$

You ask about the distribution of an order statistic $X_{(k)}$ of $N$ independent and identically distributed random variables $X_1, X_2, \ldots, X_N$ and $k=N/2$ when $N$ is even and $k=(N+1)/2$ when $N$ is odd. Because this order statistic is the median when $N$ is odd and otherwise is a median, I will refer to it henceforth as the "median."

When the common distribution function $F$ has a density $F^\prime = f,$ the order statistic has a density, too, given by

$$f_{k;X}(x) = \binom{N}{k-1;1;N-k}\, F^{k-1}(x) f(x) (1-F(x))^{N-k}.$$

By choosing time units in which the half-life is $1$, the exponential distribution function is $F(x) = 1 - e^{-x},$ whence $f(x) = e^{-x}$ and therefore

$$f_{k;X}(x) = \frac{N!}{(k-1)!(N-k)!} e^{-(N-k+1)x}\left(1 - e^{-x}\right)^{k-1}.$$

One way to understand this is to express it in terms of $U = e^{-X}$ (which, incidentally, has a Uniform$[0,1]$ distribution). The density of $U$ is

$$f_{k;U}(u) \ \propto\ u^{N-k}(1-u)^{k-1},$$

immediately recognizable as that of a Beta$(N+1-k, k)$ distribution. Thus, $X_{(k)}$ is seen as the (negative) logarithm of a Beta variable.

Here, for example, are histograms of the median of $N=12$ Exponential variables (from 40,000 simulations) and its negative exponential on which, in red, graphs of $f_{k;X}$ and $f_{k;U}$ have been overlaid to demonstrate the correctness of this analysis.

Figure

At this point you can apply what is known about the Beta distribution to determine in great detail anything you like about the distribution of $X_{(k)},$ so I won't belabor the point. Instead, consult the references, including the following CV threads:

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.