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

【问题标题】：Random numbers with user-defined continuous probability distribution具有用户定义的连续概率分布的随机数
【发布时间】：2021-03-30 16:27:33
【问题描述】：

我想模拟一些关于光子-光子-相互作用的主题。特别是，有 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 的无理数。但是，如果使用连续函数，这样的命令会更优雅。

【问题讨论】：

如果你有 CDF 的逆，你可以在[0,1) 中提取一个随机数，并使用它来获取[0, 2pi) 范围内的样本。更多信息here
您可以通过扩展 scipy.stats.rv_continuous 创建自定义分布，然后使用它来获取随机变量：docs.scipy.org/doc/scipy/reference/generated/…
这能回答你的问题吗？ How to generate random numbers with predefined probability distribution?
谢谢大家，我想我会看看逆变换采样的东西。函数 'scipy.stats.rv_continuous' 仅适用于某些概率分布。不是随意的，对吧？马尔可夫链蒙特卡罗采样也可能是某种东西。谢谢你。

标签： python random

【解决方案1】：

假设您拥有的密度函数与概率密度函数 (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

另请参阅我关于随机化的文章中的“Sampling from an Arbitrary Distribution”部分。

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

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

【讨论】：

如果我理解正确的话。你创建一个均匀分布的随机值，然后根据分布以一定的概率删除这个值。
我只是怀疑您是否以这种方式得出正确的分布函数。例如，在此分布下，您获得介于 0 和 Pi/8 之间的值（只是一个示例）的给定概率应该是 0.0788...（在 P(x) 上从 0 到 Pi/8 的积分（P(x) 是分布函数））。
x 等分布给出了一个介于 0 和 Pi/8 之间的值，概率为 0.0625。如果 x 值在此范围内，则生成的 y 值将比 x 处的计算函数值大 0.00199（在 Pmax-P(x) 上从 0 到 Pi/8 的积分）（此处 Pmax 为 64/(99 *pi))。这意味着，您的变体将输出 0 到 Pi/8 之间的值，概率为 0.0625*(1-0.00199)=0.00624 而不是 0.0788。还是我的想法有误？
即使这是正确的，也感谢您的帮助。这个想法原则上还不错。
不，它不是这样工作的。想象一个包含密度函数 (PDF) 的盒子。这个盒子是2*pi 单位宽和64/(99*pi) 单位高。 PDF 完全包含在盒子中，覆盖了盒子面积的 99/128。扔飞镖并假设它落在盒子里。那么x 和y 是飞镖的坐标。现在我们检查 dart 是否也在 PDF 中。为此，我们在x 处计算PDF 的“高度”，调用高度pdf，然后检查if y < pdf。 ...

【解决方案2】：

我有两个建议。逆变换采样方法和“删除方法”（我就这么称呼它）。逆变换采样方法：我的分布有一个反函数。但是由于域的原因，我在 math. 函数的几个地方遇到了问题。例如。 math.sqrt(-1)。您仍然需要在这里使用 if 查询。这就是为什么我决定使用 Peter 的建议。

如果您在循环中收集值并将它们绘制在直方图中，它看起来也相当不错。这里有 40000 个值和 100 个 bin

这是给感兴趣的人的完整代码

import numpy as np
import math
import random
import matplotlib.pyplot as plt

N=40000
bins=100

def Deletion_method():
    x=None
    while x==None:
        mn=0 # Lowest value of domain
        mx=2*math.pi # Highest value of domain
        bound=64/(99*math.pi) # Upper bound of PDF value

        # Choose an X inside the desired sampling domain.
        xrad=random.uniform(mn,mx)
        # Choose a Y between 0 and the maximum PDF value.
        y=random.uniform(0,bound)

        # Calculate PDF
        P=((3+(math.cos(xrad))**2)**2)*(1/(99*math.pi/4))
 
        # Does (x,y) fall in the PDF?
        if y<P:
           x=xrad
    return(x)


Values=[]

for k in range(0, N):
    Values=np.append(Values, [Deletion_method()])
   
plt.hist(Values, bins)
plt.show()

【讨论】：