How to create a random variable Rician?

I am trying to simulate a signal detection problem using Sympy and need two random variables. One with Rayleigh distribution for modeling noise and one with Rician distribution for modeling signal + noise. Sympy provides a Rayleigh distribution , but not a Rician - or at least not one of them.

What is the best way to create it? Does it exist under a different name? Is there a way to manipulate existing distributions in Rician?


Following @asmeurer's advice, I implemented my own Rice distribution, for example:

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

The distribution is similar to Wikipedia and Scipy , but it is strange that I have different results than Scipy. I will ask this question separately ( asked and answered ).

As an additional note, the following line allows you to go around the density function, which includes the Bessel function:

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

It is not generalized to all Bessel functions, but works for a modified first-order Bessel of the first type used in the Rich distribution.

+2
source share
2 answers

If you know the pdf function, it's easy to create a new distribution with sympy.stats. Take a look at existing distributions at the sympy source . You just need to subclass SingleContinuousDistribution and define some methods. For example, here is the normal distribution (with the removal of docstrings):

 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)) 
+3
source

Yes, you can create rice from a chi-square and Poisson. See a detailed discussion of rice, for example https://en.wikipedia.org/wiki/Rice_distribution :

Another case where Rice (nu, sigma) comes from the following steps:

  • We form P having a Poisson distribution with a parameter (also the average for Poisson) lambda = nu ^ 2 / (2 * sigma ^ 2).
  • Generate an X having a chi-square distribution with 2P + 2 degrees of freedom.
  • Set R = Sigma * sqrt (X).
0
source

All Articles