CODEWITHSUNDEEP
pythonscipystatistics
Denis Korzhenkov
Denis Korzhenkov

Reputation: 383

Custom distribution in scipy with pdf given

I try to define a custom distribution with pdf given via scipy.stats

import numpy as np
from scipy.stats import rv_continuous

class CustomDistribution(rv_continuous):
    def __init__(self, pdf=None):
        super(CustomDistribution, self).__init__()
        self.custom_pdf = pdf
        print "Initialized!"

    def _pdf(self, x, *args):
        if self.custom_pdf is None:
            # print 'PDF is not overridden'
            return super(CustomDistribution, self)._pdf(x, *args)
        else:
            # print 'PDF is overridden'
            return self.custom_pdf(x)

def g(x, mu):
    if x < 0:
        return 0
    else:
        return mu * np.exp(- mu * x)

my_exp_dist = CustomDistribution(pdf=lambda x: g(x, .5))
print my_exp_dist.mean()

As you see I try to define exponential distribution wuth parameter mu=0.5, but the output is as follows.

Initialized!

D:\Anaconda2\lib\site-packages\scipy\integrate\quadpack.py:357:

IntegrationWarning: The algorithm does not converge. Roundoff error is detected in the extrapolation table. It is assumed that the requested tolerance cannot be achieved, and that the returned result (if full_output = 1) is the best which can be obtained.
warnings.warn(msg, IntegrationWarning)

D:\Anaconda2\lib\site-packages\scipy\integrate\quadpack.py:357:

IntegrationWarning: The maximum number of subdivisions (50) has been achieved.

2.0576933609

If increasing the limit yields no improvement it is advised to analyze the integrand in order to determine the difficulties. If the position of a local difficulty can be determined (singularity, discontinuity) one will probably gain from splitting up the interval and calling the integrator on the subranges. Perhaps a special-purpose integrator should be used. warnings.warn(msg, IntegrationWarning)

What should I do to improve this?

NOTE: The problem of computational accuracy is discussed in this GitHub issue.

Upvotes: 3

Views: 5579

Answers (1)

Bill Bell
Bill Bell

Reputation: 21663

This seems to do what you want. An instance of the class must be given a value for the lambda parameter each time the instance is created. rv_continuous is clever enough to infer items that you do not supply but you can, of course, offer more definitions that I have here.

from scipy.stats import rv_continuous
import numpy

class Neg_exp(rv_continuous): 
    "negative exponential"
    def _pdf(self, x, lambda):
        self.lambda=lambda
        return lambda*numpy.exp(-lambda*x)
    def _cdf(self, x, lambda):
        return 1-numpy.exp(-lambda*x)
    def _stats(self,lambda):
        return [1/self.lambda,0,0,0]

neg_exp=Neg_exp(name="negative exponential",a=0)

print (neg_exp.pdf(0,.5))
print (neg_exp.pdf(5,.5))

print (neg_exp.stats(0.5))

print (neg_exp.rvs(0.5))

Upvotes: 4

Related Questions