Radioactive Decay Process Please bear with my long winded description. 

It is classical to model radioactive decay of some particles $A \rightarrow A^\ast$ by the differential equation
 \begin{align}
 \frac{\mathrm dN_A}{\mathrm dt}=-\lambda N_A
 \end{align}
where $N_A$ refers to the number of particles in state $A$ and $\lambda$ is the rate of decay. If we initially start with
$N_0$ number of particles in state $A$ then by solving the above differential equation we see that
\begin{align}
N_A(t) = N_0\mathrm e^{-\lambda t}
\end{align}
gives us an explicit number of particles in state $A$ as a function of time $t$. In particular, we also see that
\begin{align}
\lambda = \frac{\ln 2}{t_{1/2}}
\end{align}
where $t_{1/2}$ is the half-life. However, the model is completely deterministic, meaning no probability required.

However, quoting Wikipedia


  Radioactive decay is a stochastic (i.e. random) process at the level of single atoms, in that, according to quantum theory, it is impossible to predict when a particular atom will decay, regardless of how long the atom has existed. 


But collectively, the mean number of particles in state $A$ is well modeled by the above deterministic model. However, the
difficulty here is to articulate the connection between the deterministic model and the stochastic model in a rigorous manner. 

Let $T$ be a continuous positive random variable (RV) that denotes the lifetime of a single (microscopic) particle in state $A$. It is natural to make the assumption that in the stochastic model the decay of a particle in state $A$ is independent of the time that the particle spends in state $A$. With this assumption and a little bit of math (which is not important for the discussion), one can show that $T$ has an exponential distribution with a parameter $\lambda$, i.e. the
  PDF of $T$ is given by
  \begin{align}
  p_\lambda (t) =
  \begin{cases}
  \lambda\mathrm e^{-\lambda t} &amp; \text{ if } t \geq 0,\\
  0 &amp; \text{ if } t&lt;0.
  \end{cases}
  \end{align}
  and the cumulative distribution function (CDF) is given by
  \begin{align}
  P(t\geq T) = 1-\mathrm e^{-\lambda t}
  \end{align}
  for $t \geq 0$ and 0 otherwise. In particular, we see that the probability of the particle in the decayed state $A^\ast$ at time $t$ 
  is precisely $1-\mathrm e^{-\lambda t}$ and the probability in of the particle in the state $A$ is $\mathrm e^{-\lambda t}$. Hence at each fixed time
  $t$ we have a Bernoulli distribution. If we make another assumption that a collection of the radioactive particles decay independently, then we see that given $N_0$ particles at time $t$ the probability that $N$ particles survives and $N_0-N$ particles decay is given by
  the binomial distribution
  \begin{align}
  P(N_t=N) = \frac{N_0!}{N!(N_0-N)!}(1-\mathrm e^{-\lambda t})^{N_0-N}\mathrm e^{- \lambda N t}
  \end{align}
  By standard computations, we also see that the mean of the binomial distribution is $\operatorname{E}[N_t]= N_0\mathrm e^{-\lambda t}$ and the variance is $\operatorname{Var}[N_t]= N_0\mathrm e^{-\lambda t}(1-\mathrm e^{-\lambda t})$ which leads to the relative standard deviation 
  \begin{align}
 \frac{\sqrt{\operatorname{Var}[N_t]}}{\operatorname{E}[N_t]} = \frac{1}{\sqrt{N_0}}\left(\mathrm e^{\lambda t}- 1\right)^{1/2}.
  \end{align}

I have two questions.


  Question 1: Looking at the relative standard deviation, it seems to me that for very large time $t$ the deterministic model seems to be a poor approximation of the stochastic model since the relative standard deviation grows exponentially in time. Is my assertion valid? 


I actually disagree with the above calculation but have no valid argument against it. 


  Question 2: My interest in the above analysis comes from my desired to give a rigorous justification for using Monte-Carlo simulation of radioactive decay to arrive at the above deterministic model. However, I'm stuck after the above analysis. Could someone provide me with some reference or explanation to rigorous justify the connection?
The average disintegration rate is given by $\bar n=\sum_{n=0}^{N_0} nP(n)$ where $P$ is the binomial distribution. If $M$ is the average number of events in time $t$ then $\displaystyle M=N_0(1-e^{-\lambda t})$. For times smaller than the lifetime $1/\lambda$  the $\displaystyle e^{-\lambda t} \to 1-\lambda t$ then $M=N_0\lambda t$ and as the disintegration rate is $M/t$ this leads to $\displaystyle \frac{dN}{dt}=-\lambda N$. Thus provided the measuring time-bins are small compared to the lifetime this rate law should hold. 

As the decay proceeds the number of counts becomes small, perhaps $10^6$ at short times to just 1 or 2 at long times then of course the signal to noise becomes huge, zero counts in some time bins for example and 1 or 2 counts in others. (This becomes an issue in analysing data, see Anal. Chem., 2003, 75 (16), pp 4182–4187).

The expected variance is, as you quote, $V=\sigma^2=N_0e^{-\lambda t}(1-e^{-\lambda t})$ and can be rewritten as $V=Me^{-\lambda t}$ where $M$ is the mean number of events recorded in time $t$. As the time bins are normally short so that $\lambda t &lt;&lt;1$ then $\sigma=\sqrt{M}$ so that $\displaystyle \frac{\sigma}{M}=\frac{1}{\sqrt{M}}$