Overflow error subclassing a distribution using scipy.stats.rv_continuous

Question

In the documentation page of rv_continuous we can find a 'custom' gaussian being subclassed as follows.

from scipy.stats import rv_continuous
import numpy as np

class gaussian_gen(rv_continuous):
    "Gaussian distribution"
    def _pdf(self, x):
        return np.exp(-x**2 / 2.) / np.sqrt(2.0 * np.pi)
gaussian = gaussian_gen(name='gaussian')

In turn, I attempted to create a class for an exponential distribution with base 2, to model some nuclear decay:

class time_dist(rv_continuous):
    def _pdf(self, x):
        return 2**(-x)
    
random_var = time_dist(name = 'decay')

This had the purpose of then calling random_var.rvs() in order to generate a randomly distributed sample of values according to the pdf I defined. However, when I run this, I receive an OverflowError, and I don't quite understand why. Initially I thought it had to do with the fact that the function was not normalized. However, I keep making changes to the _pdf definition to no avail. Is there anything wrong with the code or is this method just ill-advised for defining functions of this sort?

Did you look into [numpy's exponential](https://numpy.org/doc/stable/reference/random/generated/numpy.random.exponential.html) and [scipy.stats.expon](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.expon.html)? You'd need `scale=1/np.log(2)` — JohanC, Dec 17 '20 at 14:32
See also [Custom distribution in scipy with pdf given](https://stackoverflow.com/questions/39494046/custom-distribution-in-scipy-with-pdf-given) — JohanC, Dec 17 '20 at 14:37

score 0 · Accepted Answer · answered Dec 17 '20 at 14:13

0

According to wikipedia, the pdf of an exponential distribution would be:

lambda * exp(-lambda*x) for x >= 0
0 for x < 0

So, probably the function should be changed as follows:

from scipy.stats import rv_continuous
import numpy as np
import matplotlib.pyplot as plt

class time_dist(rv_continuous):
    def _pdf(self, x):
        return np.log(2) * 2 ** (-x) if x >= 0 else 0

random_var = time_dist(name='decay')
plt.hist(random_var.rvs(size=500))
plt.show()

answered Dec 17 '20 at 14:13

JohanC

71,591
8
33
66

This works, thank you. Since this also inherits the `scale` and `loc` args, I can customize it further with a half life. – JTFreitas Dec 26 '20 at 21:48

Overflow error subclassing a distribution using scipy.stats.rv_continuous

1 Answers1