title: Why doesn't my Rician in Sympy match the Rician in Scipy?
body: <p>I've tried to create a Sympy continuous random variable with a Rician distribution.  With thanks to help from an <a href="https://stackoverflow.com/q/32981363/2571165">earlier question</a>, it seems that the best approach is to subclass <code>SingleContinuousDistribution</code>.  I've implemented a distribution that appears to be in agreement between <a href="https://en.wikipedia.org/wiki/Rice_distribution" rel="nofollow noreferrer">Wikipedia</a> and <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rice.html#scipy.stats.rice" rel="nofollow noreferrer">Scipy</a>, however I am not getting the same results as Scipy.</p>
  
  <p>What follows is code that implements the random variable, extracts its symbolic distribution and converts it to a Numpy representation through <code>lambdify</code> then plots my distribution against the PDF of the Scipy <code>rician</code> distribution.</p>
  
  <pre><code>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');
  </code></pre>
  
  <p>I would expect the results to match, but they don't appear to:</p>
  
  <p><a href="https://i.stack.imgur.com/LjmbR.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/LjmbR.png" alt="enter image description here"></a></p>
  
  <p>Am I doing something wrong here?</p>