时间戳

首页 > 解决方案 > 具有用户定义的连续概率分布的随机数

python - 具有用户定义的连续概率分布的随机数

问题描述

我想模拟一些关于光子相互作用的主题。特别是，有 Halpern 散射。这是关于它Halpern-Streuung的德语维基百科条目。并且差分横截面具有（3 +（cos（theta））^ 2）^ 2的角度依赖性。

我想要一个 0 到 2*Pi 之间的随机数生成器，它对应于密度函数 ((3+(cos(theta))^2)^2)*(1/(99*Pi/4) ）。因此，0、Pi 和 2*Pi 周围的值应该比 Pi/2 和 3 周围的值更频繁地出现。

我已经发现有一个关于如何随机输出具有用户定义的概率值的离散值的功能numpy.random.choice(numpy.arange(1, 7), p=[0.1, 0.05, 0.05, 0.2, 0.4, 0.2])。如果没有其他情况，我可以在紧急情况下使用它。但实际上我已经想要一个连续的概率分布。

我知道即使有这样一个 Python 命令可以输入数学分布函数，它基本上也只会产生离散的值分布，因为不能表示 1 和 0 的无理数。但是，如果使用连续函数，这样的命令会更加优雅。

标签： pythonrandom

解决方案

假设您拥有的密度函数与概率密度函数 (PDF) 成比例，您可以使用拒绝抽样方法：在框中绘制一个数字，直到该框落入密度函数内。它适用于任何具有有限域的有界密度函数，只要您知道域和边界是什么（边界是f域中的最大值）。在这种情况下，界限是64/(99*math.pi)，算法的工作原理如下：

import math
import random

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

另请参阅我关于随机化的文章中的“从任意分布中抽样”部分。

下面通过显示返回样本小于 π/8 的概率来展示该方法的正确性。为了正确起见，概率应该接近 0.0788：

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

推荐阅读