python - 为什么我的 Sympy 中的 Rician 与 Scipy 中的 Rician 不匹配？

Question

我尝试创建具有莱斯分布的 Sympy 连续随机变量。感谢 earlier question 的帮助，似乎最好的方法是子类化 SingleContinouslyDistribution。我已经实现了一个似乎在 Wikipedia 之间达成一致的发行版和 Scipy ，但是我没有得到与 Scipy 相同的结果。

接下来是实现随机变量的代码，提取其符号分布并通过 lambdify 将其转换为 Numpy 表示，然后根据 Scipy rician 的 PDF 绘制我的分布> 分布。

from sympy import *
from sympy import stats
from scipy import stats as scst
import numpy as np
import matplotlib.pyplot as plt

from sympy.stats.crv_types import rv
from sympy.stats.crv import SingleContinuousDistribution

class RicianDistribution(SingleContinuousDistribution):
    _argnames=('nu','sigma')
    @property
    def set(self): return Interval(0,oo)

    def pdf(self,x):
        nu,sigma=self.nu, self.sigma
        return (x/sigma**2)*exp(-(x**2+nu**2)/(2*sigma**2))*besseli(0,x*nu/sigma**2)

def Rician(name,nu,sigma):
    return rv(name,RicianDistribution,(nu,sigma))

#this line helps lambdify convert the sympy Bessel to a numpy Bessel
printing.lambdarepr.LambdaPrinter._print_besseli=(lambda self,expr: 'i0(%s)'%expr.argument)

x=Symbol('x')  #parameter for density function

sigma=3; pr=4

#create the symbolic Rician and numeric Rician
SpN=Rician('R',pr,sigma)  #signal plus noise
Rsci=scst.rice(pr,scale=sigma)

fx=lambdify(x,stats.density(SpN)(x),'numpy')

xs=np.linspace(0,25,1000)
plt.plot(xs,fx(xs),'b');
plt.plot(xs,Rsci.pdf(xs),'r');

我希望结果能够匹配，但它们似乎并不匹配:

我在这里做错了什么吗？

Answer 1

scipy.stats.rice中Rice分配的实现使用与 the wikipedia article 中描述的参数化略有不同的参数化.

为了使您的绘图一致，请更改此行

Rsci=scst.rice(pr,scale=sigma)

至

Rsci=scst.rice(pr/sigma, scale=sigma)

---

这里有一个更长的解释:

维基百科上显示的 PDF 是

scipy.stats.rice 文档字符串中的 PDF 为

但是，该公式没有显示 scipy 中所有连续分布都具有的 scale 参数。 (它也没有显示位置参数 loc，但我假设没有兴趣使用具有非零位置的莱斯分布。)为了创建包含尺度参数的公式，我们使用 PDF 的标准比例系列:

所以 scipy PDF 实际上是

如果我们进行识别

我们获得了维基百科文章中显示的 PDF。

在你的代码中，你的参数pr是ν，所以要转换为scipy的参数化，你必须使用b = pr/sigma。