Python：如何定义自定义发行版？

Question

我想创建一个基于 Levy 截断定律的定制分布，其内容为

p(r) = (r + r0)**(-beta)*exp(-r/k)

.

所以我用以下方式定义它：

import numpy as np
import scipy.stats as st

class LevyPDF(st.rv_continuous):
    def _pdf(self,r):
        r0 = 100
        k = 1500
        beta = 1.6
        return (r + r0)**(-beta)*np.exp(-r/k)

假设我想求

r = 0

和

r = 50km

之间的距离分布。然后：

nmin = 0
nmax = 50
my_cv = LevyPDF(a=nmin, b=nmax, name='LevyPDF')
x = np.linspace(nmin, nmax, (nmax-nmin)*2)

我不明白为什么：

sum(my_cv.cdf(x)) = 2.22

而不是

1

。

那么如何根据我定义的分布定义

N = 2000000

随机距离的直方图？

Answer 1

使用你的最小示例（稍作修改）：

import scipy.stats as st
import numpy as np
import matplotlib.pyplot as plt

class LevyPDF(st.rv_continuous):
    def _pdf(self,r):
        r0 = 100
        k = 1500
        beta = 1.6
        return (r + r0)**(-beta)*np.exp(-r/k)

nmin = 0
nmax = 50
my_cv = LevyPDF(a=nmin, b=nmax, name='LevyPDF')

要从随机变量中采样，请使用

rvs()

类中的

rv_continuous

方法：

N = 50000
X = my_cv.rvs(size=N, random_state=1)

将返回一个大小为

(N,)

的数组，其中包含从分布中采样的随机变量。使用

random_state

选项冻结您的示例并使您的脚本可重复（它定义了采样的随机种子）。

请注意，随着

N

缓慢增加，计算时间急剧增加。

要绘制直方图，请使用

matplotlib

库，请参阅

hist

:

fig, axe = plt.subplots()
n, bins, patches = axe.hist(X, 50, normed=1, facecolor='green', alpha=0.75)
plt.show(axe)

下面是从 40 自由度卡方采样的示例：

from scipy import stats
import numpy as np
import matplotlib.pyplot as plt
rv = stats.chi2(40)
N = 200000
X = rv.rvs(size=N, random_state=1)
fig, axe = plt.subplots()
n, bins, patches = axe.hist(X, 50, normed=1, facecolor='green', alpha=0.75)
plt.show(axe)

它导致：

Answer 2

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

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

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

以下是它如何解决您的问题。请注意，这比依赖 scipy.stats.rv_continuous 中的默认方法（将 pdf 与四边形积分，然后使用寻根算法反转 cdf）要快得多，后者速度较慢。我做了一个快速测试：使用 jlandercy 帖子中提出的方法生成 1000 个 rvs 需要几秒钟，但在我的计算机上使用下面的快速反演方法则不到一毫秒。

import numpy as np
from scipy.stats.sampling import NumericalInversePolynomial
from matplotlib import pyplot as plt
from scipy.integrate import quad


class MyDist:
    def __init__(self, r0, k, beta):
        self.r0 = r0
        self.k = k
        self.beta = beta

    def support(self):
        # distribution restricted to 0, 50 as indicated in the question
        # could be changed to other intervals
        return (0, 50)

    def pdf(self, x):
        # this is not a proper pdf, the normalizing
        # constant is missing (does not integrate to one)
        # this is ok for the method presented
        return (x + self.r0)**(-self.beta)*np.exp(-x/self.k)


r0 = 100
k = 1500
beta = 1.6

dist = MyDist(r0, k, beta)
gen = NumericalInversePolynomial(dist)

# compute the missing normalizing constant to plot the pdf
const_pdf = quad(dist.pdf, *dist.support())[0]

r = gen.rvs(size=10000)
x = np.linspace(r.min(), r.max(), 500)

# show histogram together with the pdf
plt.plot(x, dist.pdf(x) / const_pdf)
plt.hist(r, density=True, bins=50)
plt.show()