Drawing random samples from any distribution

Question

I want to draw random samples from any given distribution. If I follow the rv_continuous documentation, I am able to draw samples from a relatively "easy" defined function such as

class gaussian(ss.rv_continuous):

def _pdf(self, x, mu, sigma, const):
    chi2 = ((x - mu)/sigma)**2
    un_normalised = np.exp(-0.5 * chi2)
    return un_normalised/const

where const is simply the normalisation constant found using quad or simps. Now, I have

class omegaDistribution(ss.rv_continuous):

def _pdf(self, x, phi, const):
    chi2  = np.zeros(len(x))
    for i in range(len(x)):
        chi2[i] = np.dot((b - sin([x[i], phi[0]])/sigma).T, (b - sin([x[i], phi[0]]))/sigma)    
    pdf_omega = np.exp(-0.5 * chi2) * omega_prior.pdf(x)
    return pdf_omega/const

for which I have to find distribution of omega (which is denoted by x in the above code). I get the pdf as shown below:

but if I try to draw a random sample, it says:

Traceback (most recent call last):   File "non_linear.py", line 84, in <module>
    samples = omega_dist.rvs(phi = 0.02, const = norm_const, size = 1); print samples   File "/Users/Harry/anaconda/lib/python2.7/site-packages/scipy/stats/_distn_infrastructure.py", line 950, in rvs
    vals = self._rvs(*args)   File "/Users/Harry/anaconda/lib/python2.7/site-packages/scipy/stats/_distn_infrastructure.py", line 885, in _rvs
    Y = self._ppf(U, *args)   File "/Users/Harry/anaconda/lib/python2.7/site-packages/scipy/stats/_distn_infrastructure.py", line 898, in _ppf

. . . .

    chi2  = np.zeros(len(x))TypeError: object of type 'float' has no len()

Note that the distribution does not necessarily be a Gaussian. An alternative way would be to adopt a kind of accept-reject sampling method to draw random samples, but rv_continuous seems to be a much more general way to proceed.

The way I sample is as follows:

newParam = [1.0, 0.02]
norm_const = normalisation_omega(newParam, np.linspace(0.8, 1.2, 1E3))
omega_dist = omegaDistribution(name='omega_distribution', a = 1.0-0.2, b = 1.0+0.2)
x_range = np.linspace(0.8, 1.2, 1E3)
yVal =  omega_dist.pdf(x_range, phi = 0.02, const = norm_const)
samples = omega_dist.rvs(phi = 0.02, const = norm_const, size = 1E3)

Does it give an error message too, or just the traceback? That's kinda strange, if there's no message or error type. — B. Eckles, Oct 26 '16 at 13:29
Also, how are you drawing the random sample? Can you show that code too? — B. Eckles, Oct 26 '16 at 13:31
`newParam = [1.0, 0.02] # omega and phi for normalisation norm_const = normalisation_omega(newParam, np.linspace(0.8, 1.2, 1E3)) # Nomalisation constant omega_dist = omegaDistribution(name='omega_distribution', a = 1.0-0.2, b = 1.0+0.2) x_range = np.linspace(0.8, 1.2, 1E3) # Vary omega on a grid yVal = omega_dist.pdf(x_range, phi = 0.02, const = norm_const) # PDF samples = omega_dist.rvs(phi = 0.02, const = norm_const, size = 1E3) # samples` — Harry, Oct 26 '16 at 13:35
Hi Eckles...I've made a few changes to my question..would be great of you could help! Thanks. As far as I understand, the loop is the issue - but I can't find a way to get around it... — Harry, Oct 26 '16 at 13:50
The problem is that the `x` being passed into the `_pdf` function is supposed to be a float. For some reason your code assumes it's a list or other collection, and you try to run `len` on it. Is there a reason for this? — B. Eckles, Oct 26 '16 at 13:52
I've written a document about how to use rv_continuous at http://stackoverflow.com/documentation/scipy/6873/rv-continuous-for-distribution-with-parameters/23293/negative-binomial-on-positive-reals#t=20161026155219034176. If it doesn't make it clear what you have to do I'd be grateful if you would tell me so that I can improve it. — Bill Bell, Oct 26 '16 at 15:54
Hi @Bill Bell. Thanks for your reply. I tried several times, but I don't understand why it was not working in my case...I ended up generating samples using interpolation (see my blog post at [link](https://harry45.github.io/blog/2016/10/Sampling-From-Any-Distribution)). Appreciate your help! — Harry, Oct 28 '16 at 17:33
The probability density function has a discontinuity at x=0. Thus the cumulative density function is represented by an improper integral. — Bill Bell, Oct 28 '16 at 19:16

Drawing random samples from any distribution

0 Answers0