33
$\begingroup$

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.

$\endgroup$
1
  • $\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

34
$\begingroup$

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.

$\endgroup$
8
  • $\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
18
$\begingroup$

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.

$\endgroup$
1
  • 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$
    – user782
    Mar 20, 2012 at 21:39

Your Answer

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