Is there a way, using some established Python package (e.g. SciPy) to define my own probability density function (without any prior data, just $f(x) = a x + b$), so I can then make calculations with it (such as obtaining the variance of the continuous random variable)? Of course I could take, say, SymPy or Sage, create a symbolic function and do the operations, but I'm wondering whether instead of doing all this work myself I can make use of an already-implemented package.

  • $\begingroup$ Thanks for an easy way! How do you generate a histogram of random numbers implementing this way of defining custom distribution function? $\endgroup$ Jul 31, 2017 at 17:44

2 Answers 2


You have to subclass the rv_continuous class in scipy.stats

import scipy.stats as st

class my_pdf(st.rv_continuous):
    def _pdf(self,x):
        return 3*x**2  # Normalized over its range, in this case [0,1]

my_cv = my_pdf(a=0, b=1, name='my_pdf')

now my_cv is a continuous random variable with the given PDF and range [0,1]

Note that in this example my_pdf and my_cv are arbitrary names (that could have been anything), but _pdf is not arbitrary; it and _cdf are methods in st.rv_continuous one of which must be overwritten in order for the subclassing to work.

  • $\begingroup$ @GertVdE: What does "self" in def _pdf do?? $\endgroup$
    – Srivatsan
    Jun 13, 2014 at 12:40
  • $\begingroup$ @Srivatsan: see the answer on Stack Overflow $\endgroup$
    – GertVdE
    Jun 13, 2014 at 13:39
  • $\begingroup$ There is a problem with the normalization, here: you need to give a normalized probability distribution function (3*x**2, here), or the resulting random variable yields incorrect results (you can check my_cv.median(), for example). I fixed the code. $\endgroup$ Feb 23, 2016 at 17:32
  • $\begingroup$ @EOL i'm finding your use of the term "normalized" confusing. what's needed, i believe, is for the function to be centered at 0 and scaled to 1. but this answer seems to imply that the normalization needs to be over the range of x [0, 1]. can you clarify? $\endgroup$
    – abcd
    Feb 24, 2016 at 19:32
  • 1
    $\begingroup$ Maybe the standard way is to use my_cv.rvs() (which can take a size argument, for obtaining multiple samples in one go). This is what I guess from the documentation (docs.scipy.org/doc/scipy/reference/generated/…). $\endgroup$ Feb 24, 2016 at 22:05

You should check out sympy.stats. It provides an interface to deal with random variables. The following example provides a random variable X defined on the unit interval with density 2x

In [1]: from sympy.stats import *
In [2]: x = Symbol('x')
In [3]: X = ContinuousRV(x, 2*x, Interval(0, 1))

In [4]: P(X>.5) 
Out[4]: 0.750000000000000

In [5]: Var(X) # variance
Out[5]: 1/18

In [6]: E(2*cos(X)+X**2) # complex expressions are ok too
Out[6]: -7/2 + 4⋅cos(1) + 4⋅sin(1)

If you're interested this abstraction can handle some fairly complex manipulations.

  • 1
    $\begingroup$ Wow... this is just awesome! Thank you very much for this contribution. I'll keep an eye on this and your blog $\endgroup$ Mar 20, 2012 at 21:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.