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$
$\endgroup$
1
asked Mar 16, 2012 at 18:40
-
$\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 AgrawalJul 31, 2017 at 17:44
Add a comment
|
2 Answers
$\begingroup$
$\endgroup$
8
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$ Jun 13, 2014 at 12:40
-
$\begingroup$ @Srivatsan: see the answer on Stack Overflow $\endgroup$– GertVdEJun 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 checkmy_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$– abcdFeb 24, 2016 at 19:32 -
1$\begingroup$ Maybe the standard way is to use
my_cv.rvs()(which can take asizeargument, 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
$\begingroup$
$\endgroup$
1
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