Generating samples from a piecewise probability density function

586 Views Asked by At

I have the follow probability distribution function (pdf)

$$ p(x) = \begin{cases} C(x/x_{min})^{-\alpha} \quad & x \geq x_{min}\\ Ce^{-\alpha(x/x_{min}-1)} \quad & x < x_{min} \end{cases} $$

where $C$ is defined such that $p(x)$ is a pdf. The lower limit is defined by $x_0$, i.e. $p(x)=0$ for $x<x_0$, and there is no upper limit. This pdf describes a function that follows an exponential distribution until $x=x_{min}$ and a power law distribution for $x>x_{min}$. The definition of this function comes from me wanting to replicate results from this paper by Clauset et al. (p. 672).

I could resort to rejection sampling, however since I do not wish to bound my domain this seems strange to do. Also, I can, in principle, sample from both distributions separately, by using the inversion method. However, if I try to do the same for the piecewise pdf, I am stuck with computing the cumulative distribution function (cdf). Specifically, the cdf $P(x)$ has different expressions depending on if $x<x_{min}$ or $x\geq x_{min}$, since

$$ P(x) = \int_{x_0}^x dx' p(x') = \int_{x_0}^{x_{min}} dx' Ce^{-\alpha(x'/x_{min} - 1)} + \int_{x_{min}}^x dx' C(x'/x_{min})^{-\alpha} \quad \text{for } x \geq x_{min}, $$

and

$$ P(x) = \int_{x_0}^{x} dx' Ce^{-\alpha(x'/x_{min} - 1)} \quad \text{for } x < x_{min}. $$

How does one generate samples from this pdf, or any other piecewise pdf?

1

There are 1 best solutions below

1
On BEST ANSWER

Let

$$ I_1=\int_{x_0}^{x_\text{min}}\exp\left(-\alpha\left(\frac x{x_\text{min}}-1\right)\right)\mathrm dx=\frac{x_\text{min}}\alpha\left(\exp\left(-\alpha\left(\frac{x_0}{x_\text{min}}-1\right)\right)-1\right) $$

and

$$ I_2=\int_{x_\text{min}}^\infty \left(\frac x{x_\text{min}}\right)^{-\alpha}\mathrm dx=\frac{x_\text{min}}{\alpha-1}\;, $$

with $C=\frac1{I_1+I_2}$. Generate a random number $t$ uniformly in $[0,1]$. If $t\lt\frac{I_1}{I_1+I_2}$, generate an exponential variate, otherwise a power-law variate.

The standard way to generate an exponential variate over $[0,\infty]$ from a random number $r$ uniformly distributed over $[0,1]$ is $-\frac{\log r}\lambda$, where $\lambda$ is the rate parameter (in your case $\alpha/x_\text{min}$). Since your exponential distribution is restricted to a range of length $x_\text{min}-x_0$, you need to rescale $r$ to $[\exp(-\lambda(x_\text{min}-x_0)),1]=:[r_0,1]$. But you already have a random number $s$ uniformly distributed over $\left[0,\frac{I_1}{I_1+I_2}\right]$ from the first step. You can rescale that instead of generating a new one. Then your exponential variate is

$$ x_0-\frac{x_\text{min}}{\alpha}\log\left(r_0+(1-r_0)\frac{I_1+I_2}{I_1}s\right)\;. $$

The standard way to generate a power-law variate with exponent $-\alpha$ and domain $[x_\text{min},\infty]$ from a random number $r$ uniformly distributed over $[0,1]$ is $x_\text{min}r^{\frac1{1-\alpha}}$. Here, too, you already have a random number $s$, in this case uniformly distributed over $\left[\frac{I_1}{I_1+I_2},1\right]$, which you can rescale instead of generating a new one. Then your power-law variate is

$$ x_\text{min}\left(\frac{(I_1+I_2)s-I_1}{I_2}\right)^{\frac1{1-\alpha}}\;. $$