pythonmatlabnumpyrandomdistribution

equivalent to evrnd(mu,sigma,m,n) in Python


I want an alternative to this Matlab function in Python

evrnd(mu,sigma,m,n)

I think We can use something like this

numpy.random.gumbel

or just

numpy.random.uniform

Thanks in advance.


Solution

  • Matlab's evrnd generates random variates from the Gumbel distribution, also known as the Type I extreme value distribution. As explained in that link,

    The version used here is suitable for modeling minima; the mirror image of this distribution can be used to model maxima by negating R.

    You can use NumPy's implementation of the Gumbel distribution, but it uses the version of the distribution that models maxima, so you'll have to flip the values around the location (i.e. mu) parameter.

    Here's a script containing the Python function evrnd. The plot that it generates is below.

    import numpy as np
    
    
    def evrnd(mu, sigma, size=None, rng=None):
        """
        Generate random variates from the Gumbel distribution.
    
        This function draws from the same distribution as the Matlab function
    
            evrnd(mu, sigma, n)
    
        `size` may be a tuple, e.g.
    
        >>> evrnd(mu=3.5, sigma=0.2, size=(2, 5))
        array([[3.1851337 , 3.68844487, 3.0418185 , 3.49705362, 3.57224276],
               [3.32677795, 3.45116032, 3.22391284, 3.25287589, 3.32041355]])
    
        """
        if rng is None:
            rng = np.random.default_rng()
        x = -rng.gumbel(loc=-mu, scale=sigma, size=size)
        return x
    
    
    if __name__ == '__main__':
        import matplotlib.pyplot as plt
    
        mu = 10
        sigma = 2.5
        n = 20000
    
        x = evrnd(mu, sigma, n)
    
        # Plot the normalized histogram of the sample.
        plt.hist(x, bins=100, density=True, alpha=0.7)
        plt.grid(alpha=0.25)
        plt.show()
    

    plot


    If you are already using SciPy, an alternative is to use the rvs method of scipy.stats.gumbel_l. The SciPy distribution scipy.stats.gumbel_l implements the Gumbel distribution for minima, so there is no need to flip the results returned by the rvs method. For example,

    from scipy.stats import gumbel_l                                      
    
    
    mu = 10
    sigma = 2.5
    n = 20000
    
    x = gumbel_l.rvs(loc=mu, scale=sigma, size=n)