为什么我的 Sympy 中的 Rician 与 Scipy 中的 Rician 不匹配? [英] Why doesn't my Rician in Sympy match the Rician in Scipy?

问题描述

我尝试创建一个具有 Rician 分布的 Sympy 连续随机变量.感谢

我在这里做错了吗?

解决方案

Rice 发行版在

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

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

所以 scipy PDF 实际上是

如果我们进行鉴定

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

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

I've tried to create a Sympy continuous random variable with a Rician distribution. With thanks to help from an earlier question, 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, however I am not getting the same results as Scipy.

What follows is code that implements the random variable, extracts its symbolic distribution and converts it to a Numpy representation through lambdify then plots my distribution against the PDF of the Scipy rician distribution.

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');

I would expect the results to match, but they don't appear to:

Am I doing something wrong here?

解决方案

The implementation of the Rice distribution in scipy.stats.rice uses a slightly different parameterization than the parameterization described in the wikipedia article.

To make your plots agree, change this line

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

to

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

---

Here's a longer explanation:

The PDF shown on wikipedia is

The PDF in the docstring of scipy.stats.rice is

However, that formula does not show the scale parameter that all of the continuous distributions in scipy have. (It also doesn't show the location parameter loc, but I'll assume there is no interest in using a Rice distribution with a nonzero location.) To create the formula that includes the scale parameter, we use the standard scale family of a PDF:

So the scipy PDF is actually

If we make the identifications

we obtain the PDF shown in the wikipedia article.

In your code, your parameter pr is ν, so to convert to scipy's parameterization, you must use b = pr / sigma.

这篇关于为什么我的 Sympy 中的 Rician 与 Scipy 中的 Rician 不匹配?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！