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

Question

我尝试创建一个具有莱斯分布的 Sympy 连续随机变量。感谢 , it seems that the best approach is to subclass SingleContinuousDistribution. I've implemented a distribution that appears to be in agreement between Wikipedia and Scipy 的帮助，但是我没有得到与 Scipy.

相同的结果

接下来是实现随机变量的代码，提取其符号分布并通过 lambdify 将其转换为 Numpy 表示，然后根据 Scipy [=14= 的 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 uses a slightly different parameterization than the parameterization described in the wikipedia article.

中 Rice 分布的实现

为了让你的情节一致，改变这一行

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

至

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

这里有一个更长的解释：

维基百科上显示的 PDF 是

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

但是，该公式没有显示 scipy 中所有连续分布所具有的 scale 参数。（它也没有显示位置参数 loc，但我假设没有兴趣使用具有非零位置的 Rice 分布。）要创建包含比例参数的公式，我们使用标准PDF 的比例系列：

所以 scipy PDF 实际上是

如果我们进行鉴定

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

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

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

Why doesn't my Rician in Sympy match the Rician in Scipy?

python

statistics

sympy

scipy