I have a sequence of positive terms $(a_n)$, for which $\sum_{n=1}^\infty a_n = A <\infty$, and wish to take a psuedo-random sample from the discrete probability distribution
\begin{equation} \mathbf{P}[ X = n] = \frac{a_n}{A}. \end{equation}
The standard approach to sampling $X$ is to first sample a (continuous) uniform variable $U \sim \text{Unif}[0,1]$, and then set $X = x_U$, with \begin{equation} x_U = \min \Bigg\{n \, \colon \, \sum_{k=1}^n a_k > U \Bigg\}. \end{equation}
In general if the $a_n$ are arranged in decreasing order of magnitude, then for typical samples of $U$, $x_U$ can be calculated rather quickly. However, given that I want to sample the distribution many times, inevitably I will encounter cases where computing $x_U$ is very time consuming.
I am looking for an efficient implementation of this algorithm in Mathematica. My guess is that there is perhaps a way to do this with recursive functions, saving the totals of the 'large sums' (i.e. those over many $k$) to be used for future calculations.
Does anybody have any suggestions? Thanks!
EDIT
To give an example of a distirbution, consider $a_n = \log(n)/n^{(3/2)}$. This sequence is decreasing, and the associated series converges to -Zeta'[3/2]
, i.e. $A \approx 3.93224$.
On my (admittedly not particularly swift) computer, if I sample $U =0.9$, it takes approximately 1.25 seconds to return $X =2496$. My simple (naive) implementation is as follows:
X[u_] := Module[{A = -Zeta'[3/2]},
k=0;
sumk=0;
While[(sumk/Z < u),
k+=1;
sumk += N[Log[k]/k^(3/2)];
];
k];