I am trying to simulate n random discrete variable which has the following pmf

$P(X = k) = (1-p)^2kp^{k-1}$

I am thinking about using the inverse transform sampling method and I am trying to find the cdf.

$P(X \le k) = 1 - P(X > k) = 1- P(X \ge k+1) = \sum_{x=k+1}^{\infty} (1-p)^2xp^{x-1} = 1-[ (1-p)^2\sum_{x=k+1}^{\infty}xp^{x-1}]$

= $1-(1-p)^2[(k+1)p^k + (k+2)p^{k+1}+ (k+3)p^{k+2} + ...] $

What I have done so far :

Let $S = (k+1)p^k + (k+2)p^{k+1}+ (k+3)p^{k+2} + ...$

$pS =$ $(k+1)p^{k+1} + (k+2)p^{k+2}+ (k+3)p^{k+3} + ...$

then $(1-p)S = kp^k +p^k + p^{k+1} + p^{k+2}+p^{k+3}+...$ (geometric series )

$(1-p)S = kp^k + \frac{p^k}{1-p}$

$ F(x) = 1- (xp^x(1-p) + p^x)$

Then I'm trying to find the inverse :

$1- (xp^x(1-p) + p^x) < U \le 1- ((x+1)p^{x+1}(1-p) + p^{x+1})$

$xp^x(1-p) + p^x < 1-U < (x+1)p^{x+1}(1-p) + p^{x+1}$

I'm stuck here ...

Any help or hint will be appreciated !

  • 1
    $\begingroup$ I do not think there is a nice inverse for the CDF. For discrete random variable, as the CDF is just a step function, you may simply do a summation of the pmf to obtain that, and use that to generate. For example, in your case the support is $\{1, 2, 3, \ldots\}$, you generate $U \sim \text{Uniform}(0, 1)$, then check: If $U < f_X(1)$, where $f_X$ is the pmf of $X$, then assign $X$ = 1. Else check if $U < f_X(1) + f_X(2)$, then assign $X = 2$, and so on. The speed is not too bad if $f(x)$ is mainly dominated in the first few terms, and the tail rapidly converge. $\endgroup$ – BGM Dec 5 '18 at 9:12
  • $\begingroup$ I guess I have no choice $\endgroup$ – Tataria Dec 5 '18 at 9:19

There are several ways you can randomly sample from a discrete distribution, unfortunately inverting the CDF is not one of them. The plot below was generated using the Alias Method, it is particularly efficient if you know how to use binary search trees, otherwise a simple implementation of argmin, argmax will work

enter image description here

The blue bars are a histogram of $400$ samples generated with the alias method, the red points are simply

$$ P(X = x) = (1 - p)^2 x p^{x - 1} $$


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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