如何创建Rician随机变量？

Question

我正在尝试使用Sympy为信号检测问题建模，并且需要两个随机变量。一种采用瑞利分布来建模噪声，另一种采用Rician分布来建模信号+噪声。Sympy提供了Rayleigh分布，但不提供Rician分布，或者至少没有一个同名。

创建一个的最佳方法是什么？是否以其他名称存在？有没有一种方法可以将现有发行版转换为Rician？

---

遵循@asmeurer的建议，我实现了自己的Rice发行版，如下所示：

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

该分布似乎与Wikipedia和Scipy都匹配，但是奇怪的是，我得到的结果与Scipy不同。我将单独询问该问题（询问并回答）。

作为附带说明，以下几行使lambdize密度函数（包括贝塞尔函数）成为可能：

printing.lambdarepr.LambdaPrinter._print_besseli=(lambda self,expr: 'i0(%s)'%expr.argument)

它并没有推广到所有Bessel函数，但适用于Rician分布中使用的第一种零阶修改的Bessel。

Answer 1

如果您知道pdf函数，则可以使用sympy.stats创建新的发行版很容易。看一下sympy源中的现有发行版。您只需要子类化SingleContinuousDistribution并定义一些方法。例如，这是正态分布（除去了文档字符串）：

class NormalDistribution(SingleContinuousDistribution):
    _argnames = ('mean', 'std')

    @staticmethod
    def check(mean, std):
        _value_check(std > 0, "Standard deviation must be positive")

    def pdf(self, x):
        return exp(-(x - self.mean)**2 / (2*self.std**2)) / (sqrt(2*pi)*self.std)

    def sample(self):
        return random.normalvariate(self.mean, self.std)


def Normal(name, mean, std):
    return rv(name, NormalDistribution, (mean, std))