Random numbers with user-defined continuous probability distribution

Question

I would like to simulate something on the subject of photon-photon-interaction. In particular, there is Halpern scattering. Here is the German Wikipedia entry on it Halpern-Streuung. And there the differential cross section has an angular dependence of (3+(cos(theta))^2)^2.

I would like to have a generator of random numbers between 0 and 2*Pi, which corresponds to the density function ((3+(cos(theta))^2)^2)*(1/(99*Pi/4)). So the values around 0, Pi and 2*Pi should occur a little more often than the values around Pi/2 and 3.

I have already found that there is a function on how to randomly output discrete values with user-defined probability values numpy.random.choice(numpy.arange(1, 7), p=[0.1, 0.05, 0.05, 0.2, 0.4, 0.2]). I could work with that in an emergency, should there be nothing else. But actually I already want a continuous probability distribution here.

I know that even if there is such a Python command where you can enter a mathematical distribution function, it basically only produces discrete distributions of values, since no irrational numbers with 1s and 0s can be represented. But still, such a command would be more elegant with a continuous function.

Answer 1

You can use the rejection sampling approach if your density function is proportional to a probability density function (PDF): Draw a number in a box until it is within the density function's bounds. As long as you know what the domain and bound are, it works for any bounded density function with a finite domain (the bound is the maximum value of f in the domain). The bound is 64/(99*math.pi) in this case, and the algorithm is as follows:
 

import math
import random

def samples():
    minimum =0 # Lowest value of domain
    maximum =2*math.pi # Highest value of domain
    bounds =64/(99*math.pi) # Upper bound of PDF value
    while True: # Do the following until a value is returned
       # Choose an P inside the desired sampling domain.
       p=random.uniform(minimum ,maximum)
       # Choose a z between 0 and the maximum PDF value.
       z=random.uniform(0,bounds)
       # Calculate PDF
       pdf=(((3+(math.cos(p))**2)**2)*(1/(99*math.pi/4)))
       # Does (p,z) fall in the PDF?
       if z<pdf:
           # Yes, so return p
           return p
       # No, so loop

The probability that the returned sample is smaller than π/8 is shown below to demonstrate the method's accuracy. The probability should be close to 0.0788 for correctness:

print(sum(1 if samples()<math.pi/8 else 0 for _ in range(1000000))/1000000)