The differences between the two-sided geometric distribution and the "standard" geometric distribution are: 1. The samples/results can include 0 and negative numbers, and 2. the "probability of staying at 0" is calculated differently than the other integers. The Python code (below) implements Example 2.1 from the OP's [link][1]. Specifically, given the probability mass function for a two-sided geometric distribution: $$Pr[ Z = z] = \frac{1-\alpha}{1+\alpha}\alpha^{\left | z \right |}$$ - `zeroProb()` implements $\frac{1-\alpha}{1+\alpha}$ - `np.random.geometric(1-p)` implements $\alpha^{\left | z \right |}$ - `signProb()` accounts for the absolute value in $\alpha^{\left | z \right |}$ ``` #!/usr/bin/env python import random import numpy as np # The "chance of staying put at 0" is the main difference from the # standard implementation of a geometric distribution. # zeroProb() will return 1 if you "leave 0," which effectively # "turns ON" the standard implementation in twoSidedGeoDist(). def zeroProb(p): # this is the chance of "staying put at 0" if random.random() < (1.0 - p)/(1.0 + p): return 0 else: return 1 # Coin flip to determine if the result is negative or positive. # This only applies when we "leave 0." def signProb(): if random.random() < 0.5: return -1 else: return 1 def twoSidedGeoDist(p): # (1) Did we "leave 0"? [Y=1|N=0] (3) +/- return zeroProb(p) * np.random.geometric(1-p) * signProb() # (2) "Leave 0" i.e. Standard implementation alpha = 1.0/3.0 counts = {} trials = 1000000 for i in range(trials): result = twoSidedGeoDist(alpha) if result in counts: counts[result] = counts[result] + 1 else: counts[result] = 1 for x in sorted(counts.keys()): print(str(x) + "\t" + str(counts[x]) + "\t" + str(counts[x]*1.0 / trials)) ``` [1]: https://arxiv.org/pdf/0811.2841.pdf