I'm interested in implementing an algorithm for producing random samples according to some Rician distribution (sometimes known as the Ricean distribution, see).
I intend to produce these psuedo-random numbers using the Inverse Transform Sampling technique. That is, given the Cumulative Distribution Function (CDF) $C(x)$, on some PDF $P(x)$, it is possible to pick psuedo-random numbers from $P(x)$ by first generating some $U \in \text{Unif}(0,1) $, then our random number, $n$, is given by $n = C^{-1}(U)$.
My trouble is in determining $C^{-1}(x)$ for the Rice distribution. Surprisingly, I'm unable to find it online.
What I know:
The Rice distribution is characterized by some $\sigma$ and $\nu$. It's PDF is given by:
$$\frac{x}{\sigma^2}\exp\bigg(\frac{-(x^2+\nu^2)}{2\sigma^2}\bigg)I_0\bigg(\frac{x\nu}{\sigma^2}\bigg)$$
where $I_0$ is the modified Bessel function, with order $0$.
Further, the CDF of the Rician distribution is given by:
$$1 - Q_1 \bigg(\frac{\nu}{\sigma},\frac{x}{\sigma}\bigg)$$
Where $Q_1$ is the Marcum Q-function. Expanding this, I believe we have:
$$C(x) = 1 - \exp \bigg(-\frac{\frac{\nu}{\sigma}^2+\frac{x}{\sigma}^2}{2}\bigg)\sum_{k=1}^\infty \bigg(\frac{\frac{\nu}{\sigma}}{\frac{x}{\sigma}}\bigg)^k I_k(\frac{x}{\sigma} \cdot \frac{\nu}{\sigma})$$
The issue, now, is determining $C^{-1}(x)$. Unfortunately (and not surprisingly, perhaps), WolframAlpha, etc. choke up on this quickly, and I'm unsure how to start with this (dealing with $I_0$ seems particularly tricky).
If all you want to do is generate random variates with the desired distribution, then the most efficient way to do so is to generate pairs of independent normal random variates as follows: let $$X_i \sim \operatorname{Normal}(\nu, \sigma^2), \quad Y_i \sim \operatorname{Normal}(0, \sigma^2)$$ and then $R_i = \sqrt{X_i^2 + Y_i^2}$ will obey $$R_i \sim \operatorname{Rice}(|\nu|, \sigma).$$ This is because this process is how the Rice distribution arises. See https://en.wikipedia.org/wiki/Rice_distribution#Related_distributions.
In the same section, there is another method that is described: let $$P_i \sim \operatorname{Poisson}\left(\lambda = \frac{\nu^2}{2\sigma^2}\right), \\ X_i \mid P_i \sim \chi^2 (2P_i + 2).$$ Then $R_i = \sigma \sqrt{X_i}$ is also Rice-distributed.
Either of these methods will be much faster than the inverse transform sampling, but between the two, I prefer the first. You can modify the Box-Muller method in a very simple way: just generate the pair $$(Z_0, Z_1) = \left(\sqrt{-2 \log U_1} \cos (2\pi U_2), \sqrt{-2 \log U_2} \sin (2\pi U_2)\right)$$ where $U_1, U_2$ are independent uniform on $(0,1)$, and then $$(X_i, Y_i) = (\nu + \sigma Z_0, \sigma Z_1)$$ or more directly $$\begin{align} R_i &= \sqrt{(\nu + \sigma Z_0)^2 + (\sigma Z_1)^2} \\ &= \sqrt{\nu^2 + 2 \sigma \nu \sqrt{-2 \log U_1} \cos (2 \pi U_2) - 2\sigma^2 \log U_1} \end{align}$$ will generate the desired variate.
I should report that I was able to use the first method (two normal variates via Box-Muller) successfully, but I have not been able to get the second method (chi-squared with Poisson hierarchical model) to work.