Random numbers with user-defined continuous probability distribution

0 votes

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.

Mar 23 in Machine Learning by Dev
• 6,000 points
152 views

1 answer to this question.

0 votes

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)
answered Mar 25 by Nandini
• 5,480 points

Related Questions In Machine Learning

0 votes
0 answers

Create random numbers with left skewed probability distribution

I would like to pick a number ...READ MORE

Mar 25 in Machine Learning by Dev
• 6,000 points
374 views
+1 vote
2 answers

ValueError: Found input variables with inconsistent numbers of samples: [1, 1000]

Hi@akhtar, Here you used x as your feature ...READ MORE

answered Apr 14, 2020 in Machine Learning by MD
• 95,380 points

edited Aug 11, 2021 by Soumya 50,356 views
0 votes
1 answer

Found input variables with inconsistent numbers of samples:

Hi@sagar, You have converted your Dataframe into an ...READ MORE

answered Jul 13, 2020 in Machine Learning by MD
• 95,380 points
4,244 views
0 votes
1 answer

problem with Found input variables with inconsistent numbers of samples: [1204, 134]

Hi@Alessandro, Here you used x as your feature parameter ...READ MORE

answered Jul 20, 2020 in Machine Learning by MD
• 95,380 points
3,599 views
0 votes
2 answers
+1 vote
2 answers

how can i count the items in a list?

Syntax :            list. count(value) Code: colors = ['red', 'green', ...READ MORE

answered Jul 7, 2019 in Python by Neha
• 330 points

edited Jul 8, 2019 by Kalgi 3,128 views
0 votes
1 answer
0 votes
1 answer
0 votes
1 answer

Plotting joint probability of two random variable choices

Create pandas.DataFrame , you can use seaborn.jointplot ...READ MORE

answered Mar 23 in Machine Learning by Nandini
• 5,480 points
349 views
webinar REGISTER FOR FREE WEBINAR X
REGISTER NOW
webinar_success Thank you for registering Join Edureka Meetup community for 100+ Free Webinars each month JOIN MEETUP GROUP