11

I'm trying to create a distribution based on some data I have, then draw randomly from that distribution. Here's what I have:

from scipy import stats import numpy def getDistribution(data): kernel = stats.gaussian_kde(data) class rv(stats.rv_continuous): def _cdf(self, x): return kernel.integrate_box_1d(-numpy.Inf, x) return rv() if __name__ == "__main__": # pretend this is real data data = numpy.concatenate((numpy.random.normal(2,5,100), numpy.random.normal(25,5,100))) d = getDistribution(data) print d.rvs(size=100) # this usually fails

I think this is doing what I want it to, but I frequently get an error (see below) when I try to do d.rvs(), and d.rvs(100) never works. Am I doing something wrong? Is there an easier or better way to do this? If it's a bug in scipy, is there some way to get around it?

Finally, is there more documentation on creating custom distributions somewhere? The best I've found is the scipy.stats.rv_continuous documentation, which is pretty spartan and contains no useful examples.

The traceback:

Traceback (most recent call last): File "testDistributions.py", line 19, in print d.rvs(size=100) File "/usr/local/lib/python2.6/dist-packages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/stats/distributions.py", line 696, in rvs vals = self._rvs(*args) File "/usr/local/lib/python2.6/dist-packages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/stats/distributions.py", line 1193, in _rvs Y = self._ppf(U,*args) File "/usr/local/lib/python2.6/dist-packages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/stats/distributions.py", line 1212, in _ppf return self.vecfunc(q,*args) File "/usr/local/lib/python2.6/dist-packages/numpy-1.6.1-py2.6-linux-x86_64.egg/numpy/lib/function_base.py", line 1862, in call theout = self.thefunc(*newargs) File "/usr/local/lib/python2.6/dist-packages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/stats/distributions.py", line 1158, in _ppf_single_call return optimize.brentq(self._ppf_to_solve, self.xa, self.xb, args=(q,)+args, xtol=self.xtol) File "/usr/local/lib/python2.6/dist-packages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/optimize/zeros.py", line 366, in brentq r = _zeros._brentq(f,a,b,xtol,maxiter,args,full_output,disp) ValueError: f(a) and f(b) must have different signs

Edit

For those curious, following the advice in the answer below, here's code that works:

from scipy import stats import numpy def getDistribution(data): kernel = stats.gaussian_kde(data) class rv(stats.rv_continuous): def _rvs(self, *x, **y): # don't ask me why it's using self._size # nor why I have to cast to int return kernel.resample(int(self._size)) def _cdf(self, x): return kernel.integrate_box_1d(-numpy.Inf, x) def _pdf(self, x): return kernel.evaluate(x) return rv(name='kdedist', xa=-200, xb=200)
Improve this question
2
  • So when we are doing the above and calling randoms = getDistribution(Mydata) and then randoms = randoms.rvs(size=1000) does it perform the three def inside the class? i.e calculating pdf, integrating it etc? Commented Jul 11, 2014 at 11:07
  • I do get my randoms to follow the data distribution, but I would like to smooth it so that it does not follow the data distribution exactly. I have been manually adjusting bandwidth in kernel to do that. For example, something like how we specify a PDF function and then use the PDF function to create randoms using Metropolis Hastings. Commented Jul 11, 2014 at 11:10

1 Answer 1

Reset to default
7

Specifically to your traceback:

rvs uses the inverse of the cdf, ppf, to create random numbers. Since you are not specifying ppf, it is calculated by a rootfinding algorithm, brentq. brentq uses lower and upper bounds on where it should search for the value at with the function is zero (find x such that cdf(x)=q, q is quantile).

The default for the limits, xa and xb, are too small in your example. The following works for me with scipy 0.9.0, xa, xb can be set when creating the function instance

def getDistribution(data): kernel = stats.gaussian_kde(data) class rv(stats.rv_continuous): def _cdf(self, x): return kernel.integrate_box_1d(-numpy.Inf, x) return rv(name='kdedist', xa=-200, xb=200)

There is currently a pull request for scipy to improve this, so in the next release xa and xb will be expanded automatically to avoid the f(a) and f(b) must have different signs exception.

There is not much documentation on this, the easiest is to follow some examples (and ask on the mailing list).

edit: addition

pdf: Since you have the density function also given by gaussian_kde, I would add the _pdf method, which will make some calculations more efficient.

edit2: addition

rvs: If you are interested in generating random numbers, then gaussian_kde has a resample method. Random Samples can be generated by sampling from the data and adding gaussian noise. So, this will be faster than the generic rvs using the ppf method. I would write a ._rvs method that just calls gaussian_kde's resample method.

precomputing ppf: I don't know of any general way to precompute the ppf. However, the way I thought of doing it (but never tried so far) is to precompute the ppf at many points and then use linear interpolation to approximate the ppf function.

edit3: about _rvs to answer Srivatsan's question in the comment

_rvs is the distribution specific method that is called by the public method rvs. rvs is a generic method that does some argument checking, adds location and scale, and sets the attribute self._size which is the size of the requested array of random variables, and then calls the distribution specific method ._rvs or it's generic counterpart. The extra arguments in ._rvs are shape parameters, but since there are none in this case, *x and **y are redundant and unused.

I don't know how well the size or shape of the .rvs method works in the multivariate case. These distributions are designed for univariate distributions, and might not fully work for the multivariate case, or might need some reshapes.

Improve this answer
Sign up to request clarification or add additional context in comments.

3 Comments

Awesome, thanks, this was very helpful. Is there some way to precompute the ppf from the cdf using the same methods that scipy is using, so that it's more efficient? I notice that _cdf() gets called a lot for each rv() call.
I added a few more comments on rvs and ppf. One more comment: gaussian_kde will not be very good in the tails if you have data with fat tails. When I thought about writing a similar distribution subclass, I would have tried to use pareto tails. I read a comment about this in a forum, and matlab has a Pareto Tail Distribution.
@user333700 I would like to know what does def _rvs(self, *x, **y) in the answer posted as the edit in the question do??

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.