vendredi 3 avril 2020

How to create a random variable with a sharp custom p.d.f using scipy.stats.rv_continuous?

Context

I'm trying to use scipy.stats.rv_continuous to create a random variable with a "sharp" p.d.f: f(x) = 2*x if x is in [0,1], otherwise f(x) = 0. It throws an error which I understand to be related to the integral of this p.d.f evaluated using scipy.integrate.quad is not 1 (related to my question detailed in this thread).

I've since figured out that to fix the integral, I can simply specify the points argument to the scipy.integrate.quad function:

points(sequence of floats,ints), optional

A sequence of break points in the bounded integration interval where local difficulties of the integrand may occur (e.g., singularities, discontinuities). The sequence does not have to be sorted. Note that this option cannot be used in conjunction with weight.

However, I haven't found a way to specify this type of break points to the rv_continuous class. Hence the error remains. Any idea how to fix it? Code snippet and error message are listed below.

Code to reproduce the error

# use scipy to create random number for f(x) = 2x when x in [0,1] and 0, otherwise
from scipy.stats import rv_continuous
class custom_rv(rv_continuous):
    "custom distribution"
    def _pdf(self, x):
        if x >= 0.0 and x <=1.0:
            return 2*x
        else:
            return 0.0
rv = custom_rv(name='2x')
rvs = rv.rvs(size=1)

Output/Error

RuntimeError: Failed to converge after 100 iterations.

Code for the integral fix

from scipy.integrate import quad
print(quad(rv._pdf, -10.0, 10.0))
print(quad(rv._pdf, -10.0, 10.0, points = [1]))
# this following line throws an error though
# since it doesn't support infinity inputs combined with break points.
print(quad(rv._pdf, -np.inf, np.inf, points = [1]))

# Output:
(0.0, 0.0)
(1.0000000000149878, 2.3295111147512406e-09)
ValueError: Infinity inputs cannot be used with break points.



Aucun commentaire:

Enregistrer un commentaire