I am interested in a particular density, and I need to sample it "regularly" in a way that represent its shape (not random).
Formally, f is my density function, F is the corresponding cumulative density function (F' = f), whose reverse function rF = F^-1 does exist. I am interested in casting a regular sample from [0, 1] into my variable domain through F^-1. Something like:
import numpy as np
uniform_sample = np.linspace(0., 1., 256 + 2)[1:-1] # source sample
shaped_sample = rF(uniform_sample) # this is what I want to get
Is there a dedicated way to do this with numpy, or should I do this by hand? Here is the 'by hand' way for exponential law:
l = 5. # exponential parameter
# f = lambda x: l * np.exp(-l * x) # density function, not used
# F = lambda x: 1 - np.exp(-l * x) # cumulative density function, not used either
rF = lambda y: np.log(1. / (1. - y)) / l # reverse `F^-1` function
# What I need is:
shaped_sample = rF(uniform_sample)
I know that, in theory, rF is internally used for drawing random samples when np.random.exponential is called, for example (a uniform, random sample from [0, 1] is transformed by rF to get the actual result). So my guess is that numpy.random does know the rF function for each distribution it offers.
How do I access it? Does numpy provide functions like:
np.random.<any_numpy_distribution>.rF
or
np.random.get_reverse_F(<any_custom_density_function>)
.. or should I derive / approximate them myself?
scipy has probability distribution objects for all (I think) of the probability distributions in numpy.random.
http://docs.scipy.org/doc/scipy/reference/stats.html
The all have a ppf() method that does what you want.
http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.ppf.html
In your example:
import scipy.stats as st
l = 5. # exponential parameter
dist = st.expon(0., l) # distribution object provided by scipy
f = dist.pdf # probability density function
F = dist.cdf # cumulative density function
rF = dist.ppf # percent point function : reverse `F^-1` function
shaped_sample = rF(uniform_sample)
# and much more!