pythonmathscipystatistics

How to sample from a distribution given the CDF in Python


I would like to draw samples from a probability distribution with CDF 1 - e^(-x^2).

Is there a method in python/scipy/etc. to enable you to sample from a probability distribution given only its CDF?


Solution

  • To create a custom random variable class given a CDF you could subclass scipy.rv_continuous and override rv_continuous._cdf. This will then automatically generate the corresponding PDF and other statistical information about your distribution, e.g.

    import matplotlib.pyplot as plt
    import numpy as np
    from scipy import stats
    
    class MyRandomVariableClass(stats.rv_continuous):
        def __init__(self, xtol=1e-14, seed=None):
            super().__init__(a=0, xtol=xtol, seed=seed)
    
        def _cdf(self, x):
            return 1-np.exp(-x**2)
    
    
    if __name__ == "__main__":
        my_rv = MyRandomVariableClass()
    
        # sample distribution
        samples = my_rv.rvs(size = 1000)
    
        # plot histogram of samples
        fig, ax1 = plt.subplots()
        ax1.hist(list(samples), bins=50)
    
        # plot PDF and CDF of distribution
        pts = np.linspace(0, 5)
        ax2 = ax1.twinx()
        ax2.set_ylim(0,1.1)
        ax2.plot(pts, my_rv.pdf(pts), color='red')
        ax2.plot(pts, my_rv.cdf(pts), color='orange')
    
        fig.tight_layout()
        plt.show()