Error when drawing random numbers from a custom continuous distribution using scipy.rv_continuous

Question

Error when drawing random numbers from a custom continuous distribution using scipy.rv_continuous

pelegs

2022年4月14日 14:34

I am trying to generate a sample of random numbers from a custom distribution

$$ p(x) = x^{n}e^{-xtn}. $$

After reading the tutorial on scipy's website, I wrote a subclass which I called kbayes:

class kbayes(rv_continuous):
    def _pdf(self, x, t, n):
        p = x**n * np.exp(-t*n*x)
        s = np.sum(p)
        return p/s

The line s=np.sum(p) is there to normalize the distribution.

The pdf seems to be ok when I check it on some numbers: running the following code

ks = np.logspace(-5, -2.3, 1000)
p = kb._pdf(ks, 500, 12)
hist, _ = np.histogram(p, ks)
plt.plot(ks, p)

results in

which is what I expect.

But when I try to generate a sample, via e.g.

p_test = kb.rvs(size=1000, t=500, n=12)

I get the following error:

/home/[redacted]/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in double_scalars
  """
/home/[redacted]/anaconda3/lib/python3.6/site-packages/scipy/integrate/quadpack.py:385: IntegrationWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  warnings.warn(msg, IntegrationWarning)

I understand that this distribution has a vert small practical range (where the probability is higher than $10^{-10}$), but I expected the rv_continuous class to be able to cope with that.

I would appreciate any input on how to make this work, or any other way to draw a sample from this distribution using python.

Topic numpy scipy python statistics

Category Data Science

Brian Spiering · Accepted Answer · 2022年4月14日 14:34

That is a warning which can be ignored. Below is a cleaned-up example that generates samples and ignores that warning.

from   warnings import filterwarnings

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

filterwarnings('ignore')

class kbayes_gen(rv_continuous):
    def _pdf(self, x, t, n):
        p = x**n * np.exp(-t*n*x)
        s = np.sum(p)
        return p/s

kb = kbayes_gen()    
ks = np.logspace(-5, -2.3, 1000)
p = kb._pdf(x=ks, t=500, n=12)
plt.plot(ks, p);

samples = kb.rvs(size=1000, t=500, n=12)

Error when drawing random numbers from a custom continuous distribution using scipy.rv_continuous

About