Search code examples
pythonscipy

SciPy - Custom probability density function and generating samples from it

I have the following probability density function: f(x) = 0.5*exp(-0.5*(x - a)) where a is a parameter. For each value of a we have one distinct pdf. Following this Q&A, my code is as below:

class exp_pdf(rv_continuous):
    def _get_support(self, param):
        return 0, float('inf')
    def _pdf(self, x, param):
        return beta*np.exp(-beta*(x - param))

Then I plot the PDF for an arbitrary value of param, for example, param = 1. The plot is correct:

param_temp = 1
exp_dist = exp_pdf(a = param_temp, name='Custom Exponential Function')
points = np.linspace(1, 10, 100)
plt.plot(points, exp_dist.pdf(points, param = param_temp), 'r-', lw=2)

enter image description here

However, when I try to generate samples from this pdf, the samples are not matched:

samples = exp_dist.rvs(size = 100, param = param_temp)
plt.hist(samples, bins = 30, density = True, alpha = 0.7)

enter image description here

If I don't use _get_support and hardcode the value of param, the histogram and the density function match well:

class exp_pdf(rv_continuous):
    def _pdf(self, x):
        return beta*np.exp(-beta*(x - 1))

exp_dist = exp_pdf(a = 1, name = 'Custom Exponential Function')
points = np.linspace(1, 10, 100)
plt.plot(points, exp_dist.pdf(points), 'r-', lw = 2)

samples = exp_dist.rvs(size = 100)
plt.hist(samples, bins = 30, density = True, alpha = 0.7)

enter image description here

May I ask for your suggestions on this matter? How do I get _get_support to work correctly for this case?

Solution

  • Your distribution is not defined for x < param. Try this definition:

    class exp_pdf(rv_continuous):
    def _get_support(self, param):
        return param, float('inf')
    def _pdf(self, x, param):
        return beta*np.exp(-beta*(x - param))