How to simulate from an (arbitrary) continuous probability distribution?

Question

I have a probability density function like this:

def p1(x):
    return ( sin(x) ** (-0.75) ) / (4.32141 * (x ** (1/5)))

I want to denerate random value on [0; 1] with this pdf. How can I do random value?

Sampling from an arbitrary distribution can be done by sampling uniformly in [0, 1] then using the inverse of the cumulative density function. If you cannot compute this inverse-cdf analytically, you can use numerical integration of your pdf (e.g. using trapezoidal rule) and store the values in a list; then you can do a binary search of your [0, 1] sample to find a sample of your pdf. — Francis Colas, Mar 17 '15 at 09:22
@FrancisColas, it is a typical way to solve my task. But is it really, that I can solve my problem using classic way, when there is so powerful tool as scipy? — Denis, Mar 17 '15 at 09:26
Well, I don't know of any specific way using numpy/scipy. But notice that it is just a comment to give you a hint for a solution, others might post more satisfying answers. — Francis Colas, Mar 17 '15 at 09:41
@FrancisColas Better integrate it with `quad` then interpolate the inverse function. — ev-br, Mar 17 '15 at 10:07

score 15 · Accepted Answer · answered Mar 17 '15 at 09:59

15

As mentioned by Francis you'd better know the cdf of your distribution. Anyway scipy provides a handy way to define custom distributions. It looks pretty much like that

from scipy import stats
class your_distribution(stats.rv_continuous):
    def _pdf(self, x):
        return ( sin(x) ** (-0.75) ) / (4.32141 * (x ** (1/5)))

distribution = your_distribution()
distribution.rvs()

answered Mar 17 '15 at 09:59

Gioelelm

2,645
5
30
49

2

(+1), good answer. And if this turns out to be too slow (since it goes through the generic code paths), only then it's worth to integrate/interpolate the cdf and it's inverse for random sampling. – ev-br Mar 17 '15 at 10:06
pdf_z is an interpolation object of scipy.interpolate.interp1d, bounds=0 if outside x the range ! from scipy import stats class TransfromPriorPDF(stats.rv_continuous): def _pdf(self, x): return pdf_z(x) dist = TransfromPriorPDF() dist.rvs() This gives me a runtime error: --> 776 r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp) 777 return results_c(full_output, r) 778 RuntimeError: Failed to converge after 100 iterations. Any idea where it is going wrong? – Sebastiano1991 Feb 02 '21 at 11:01
Is this solution supportive to the N-dimensional distribution? For example, the PDF may be **exp(-1/(x^2+y^2))**, **exp(-1/(a^2+b^2+c^2+d^2))** or something like them. (The examples may be not guaranteed to be real PDFs.) I may want to sample a vector **(a,b,c,d)**. – BinChen Mar 26 '23 at 15:10

score 4 · Answer 2 · answered Mar 01 '17 at 15:16

Without using scipy and given a numerical sampling of your PDF, you can sample using a cumulative distribution and linear interpolation. The code below assumes equal spacing in x. It could be modified to do an integration for an arbitrarily sampled PDF. Note it renormalises the PDF to 1 within the range of x.

import numpy as np

def randdist(x, pdf, nvals):
    """Produce nvals random samples from pdf(x), assuming constant spacing in x."""

    # get cumulative distribution from 0 to 1
    cumpdf = np.cumsum(pdf)
    cumpdf *= 1/cumpdf[-1]

    # input random values
    randv = np.random.uniform(size=nvals)

    # find where random values would go
    idx1 = np.searchsorted(cumpdf, randv)
    # get previous value, avoiding division by zero below
    idx0 = np.where(idx1==0, 0, idx1-1)
    idx1[idx0==0] = 1

    # do linear interpolation in x
    frac1 = (randv - cumpdf[idx0]) / (cumpdf[idx1] - cumpdf[idx0])
    randdist = x[idx0]*(1-frac1) + x[idx1]*frac1

    return randdist

Good answer - I used this for an very complicated pdf (Ratio of two normals) and it works a charm (plus it is super fast) - you should post your answer in the linked SO post since this is marked as duplicate - I think your answer is easier to implement - here is the pdf: http://static.stevereads.com/papers_to_read/on_the_ratio_of_two_correlated_normal_random_variables.pdf — Xavier Bourret Sicotte, May 31 '19 at 23:15

How to simulate from an (arbitrary) continuous probability distribution?

2 Answers2