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, 2022 1,282 views

## 1 answer to this question.

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)`
• 5,480 points

## Create random numbers with left skewed probability distribution

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

+1 vote

## how do i change string to a list?

suppose you have a string with a ...READ MORE

## how can i randomly select items from a list?

You can also use the random library's ...READ MORE

+1 vote

## how can i count the items in a list?

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