20
$\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$
  • $\begingroup$ Thanks for an easy way! How do you generate a histogram of random numbers implementing this way of defining custom distribution function? $\endgroup$ – Ankur Agrawal Jul 31 '17 at 17:44
23
$\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$
  • $\begingroup$ @GertVdE: What does "self" in def _pdf do?? $\endgroup$ – Srivatsan Jun 13 '14 at 12:40
  • $\begingroup$ @Srivatsan: see the answer on Stack Overflow $\endgroup$ – GertVdE Jun 13 '14 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$ – Eric O Lebigot Feb 23 '16 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$ – dbliss Feb 24 '16 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$ – Eric O Lebigot Feb 24 '16 at 22:05
15
$\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$
  • $\begingroup$ Wow... this is just awesome! Thank you very much for this contribution. I'll keep an eye on this and your blog $\endgroup$ – astrojuanlu Mar 20 '12 at 21:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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