通过python生成具有一定概率密度函数的随机变量

Question

我想生成具有麦克斯韦能量（例如 f(ene)）或任意分布的随机变量。 概率密度函数如下：

def f(ene):  
     le=3  
     return 2(ene/pi/le**3)*np.exp(-ene/le)

我想生成 10000 个样本，例如

f.rvs(scale,size)

生成一个包含 10000 个元素的列表，其密度函数与 f(ene) 匹配

我该怎么办？

示例： http://docs.scipy.org/doc/scipy/reference/ generated/scipy.stats.maxwell.html scipy 提供麦克斯韦分布： 使用方法：rvs(loc=0,scale=1,size=1)随机变量。

print maxwell.rvs(1,10000)

将生成 10000 个具有麦克斯韦分布的样本。但scipy只提供了一些功能。现在我有另一个函数不在 scipy's 中。我怎么能这么做呢？

Answer 1

好的，让我们从分发开始吧。

速度的概率密度函数是

PDF(v) = sqrt(2/PI)*v^2*exp(-v^2/2a^2)/a^3

让我们将其从速度转换为能量，考虑 E=v^2/2（并假设质量等于 1）

PDF(E) = sqrt(2/PI)*2*E*exp(-E/a^2)/a^3

看起来还好吗？不，这是错误的。为什么？因为 PDF 乘以间隔给出了概率，而我们没有转换间隔。速度的概率为

PDF(v) * dv

，我们必须插入雅可比行列式

PDF(E)*dE = sqrt(2/PI)*2*E*exp(-E/a^2)/a^3 * |dv/dE| * dE

|dv/dE| = 1/|dE/dv| = 1/v = 1 / sqrt(2*E)

大家一起

PDF(E) = 2*sqrt(1/PI)*sqrt(E)*exp(-E/a^2)/a^3

哪里

a=sqrt(k*T)

您的发行版丢失了

sqrt(E)

我将如何采样它。我会利用它实际上是三个独立的一维分布的乘积这一事实

PDF(vx) = 1/sqrt(2*PI)*exp(-vx^2/2a^2)
PDF(vy) = 1/sqrt(2*PI)*exp(-vy^2/2a^2)
PDF(vz) = 1/sqrt(2*PI)*exp(-vz^2/2a^2)

这意味着它们都是高斯分布的。如此简单的算法：从均值等于 0 且西格玛等于

vx

的高斯样本

a

，从同一分布独立地采样

vy

，对

vz

进行样本并将它们组合在一起

v = sqrt(vx^2 + vy^2 + vz^2)

或者如果您需要采样能量

E = ( vx^2 + vy^2 + vz^2 )/ 2

Answer 2

这是我能想到的最好的办法，但我对你在这里尝试做的事情没有信心：

import numpy as np

def f(ene):
    le = 3
    return 2 * (ene / np.pi / le ** 3) * np.exp(-ene / le)

def rvs(scale, size):
    samples = np.random.sample((size,))
    return scale * f(samples)

我得到了这个结果。显然，每次都会不同，因为这是随机的。

In:  print(rvs(100, 10))
Out: [ 1.56035154  0.85509302  0.76543496  1.36966862  0.20596924  0.0395071
  1.35029318  0.69437599  0.59772671  0.56721737]

Answer 3

在过去的几年里，SciPy 中添加了一些不错的新工具来解决 Python 中的此类问题。 只需提供一些有关分布的信息（例如密度/pdf），您就可以轻松地从自定义连续或离散单变量分布生成样本。

不同方法的概述： https://docs.scipy.org/doc/scipy/reference/stats.sampling.html

教程： https://docs.scipy.org/doc/scipy/tutorial/stats/sampling.html