title: Custom numpy (or scipy?) probability distribution for random number generation
body: <h2>The issue</h2>
  <p><strong>Tl;dr</strong>: I would like a function that randomly returns a float (or optionally an ndarray of floats) in an interval following a probability distribution that resembles the sum of a &quot;Gaussian&quot; and a uniform distributions.</p>
  <p>The function (or class) - let's say <code>custom_distr()</code> - should have as inputs (with default values already given):</p>
  <ul>
  <li>the lower and upper bounds of the interval: <code>low=0.0</code>, <code>high=1.0</code></li>
  <li>the mean and standard deviation parameters of the &quot;Gaussian&quot;: <code>loc=0.5</code>, <code>scale=0.02</code></li>
  <li>the size of the output: <code>size=None</code></li>
  <li><code>size</code> can be an integer or a tuple of integers. If so, then <code>loc</code> and <code>scale</code> can either both simultaneously be scalars, or ndarrays whose <code>shape</code> corresponds to <code>size</code>.</li>
  </ul>
  <p>The output is a scalar or an ndarray, depending on size.</p>
  <p>The output has to be scaled to certify that the cumulative distribution is equal to 1 (I'm uncertain how to do this).</p>
  <p>Note that I'm following <code>numpy.random.Generator</code>'s naming convention from <code>uniform</code> and <code>normal</code> distributions as reference, but the nomenclature and the utilized packages does not really matter to me.</p>
  <h2>What I've tried</h2>
  <p>Since I couldn't find a way to &quot;add&quot; <code>numpy.random.Generator</code>'s uniform and Gaussian distributions directly, I've tried using <code>scipy.stats.rv_continuous</code> subclassing, but I'm stuck at how to define the <code>_rvs</code> method, or the <code>_ppf</code> method to make it fast.</p>
  <p>From what I've understood of <a href="https://github.com/scipy/scipy/blob/v1.6.3/scipy/stats/_continuous_distns.py" rel="nofollow noreferrer">rv_continuous class definition in Github</a>, <code>_rvs</code> uses <code>numpy</code>'s <code>random.RandomState</code> (which is out of date in comparison to <code>random.Generator</code>) to make the distributions. This seems to defeat the purpose of using <code>scipy.stats.rv_continuous</code> subclassing.</p>
  <p>Another option would be to define <code>_ppf</code>, the percent-point function of my custom distribution, since according to <a href="https://github.com/scipy/scipy/blob/v1.6.3/scipy/stats/_distn_infrastructure.py#L1458-L2627" rel="nofollow noreferrer">rv_generic class definition in Github</a>, the default function <code>_rvs</code> uses <code>_ppf</code>. But I'm having trouble defining this function by hand.</p>
  <p>Following, there is a MWE, tested using <code>low=0.0</code>, <code>high=1.0</code>, <code>loc=0.3</code> and <code>scale=0.02</code>. The names are different than the &quot;<strong>The issue</strong>&quot; section, because terminologies of terms are different between <code>numpy</code> and <code>scipy</code>.</p>
  <pre class="lang-py prettyprint-override"><code>import numpy as np
  from scipy.stats import rv_continuous
  import scipy.special as sc
  import matplotlib.pyplot as plt
  import time
  
  
  # The class definition
  class custom_distr(rv_continuous):
      def __init__(self, my_loc=0.5, my_scale=0.5, a=0.0, b=1.0, *args, **kwargs):
          super(custom_distr, self).__init__(a, b, *args, **kwargs)
          self.a = a
          self.b = b
          self.my_loc = my_loc
          self.my_scale = my_scale
  
      def _pdf(self, x):
          # uniform distribution
          aux = 1/(self.b-self.a)
          # gaussian distribution
          aux += 1/np.sqrt(2*np.pi*self.my_scale**2) * \
                   np.exp(-(x-self.my_loc)**2/2/self.my_scale**2)
          return aux/2  # divide by 2?
  
      def _cdf(self, x):
          # uniform distribution
          aux = (x-self.a)/(self.b-self.a)
          # gaussian distribution
          aux += 0.5*(1+sc.erf((x-self.my_loc)/(self.my_scale*np.sqrt(2))))
          return aux/2  # divide by 2?
  
  
  # Testing the class
  if __name__ == &quot;__main__&quot;:
      my_cust_distr = custom_distr(name=&quot;my_dist&quot;, my_loc=0.3, my_scale=0.02)
  
      x = np.linspace(0.0, 1.0, 10000)
  
      start_t = time.time()
      the_pdf = my_cust_distr.pdf(x)
      print(&quot;PDF calc time: {:4.4f}&quot;.format(time.time()-start_t))
      plt.plot(x, the_pdf, label='pdf')
  
      start_t = time.time()
      the_cdf = my_cust_distr.cdf(x)
      print(&quot;CDF calc time: {:4.4f}&quot;.format(time.time()-start_t))
      plt.plot(x, the_cdf, 'r', alpha=0.8, label='cdf')
  
      # Get 10000 random values according to the custom distribution
      start_t = time.time()
      r = my_cust_distr.rvs(size=10000)
      print(&quot;RVS calc time: {:4.4f}&quot;.format(time.time()-start_t))
  
      plt.hist(r, density=True, histtype='stepfilled', alpha=0.3, bins=40)
  
      plt.ylim([0.0, the_pdf.max()])
      plt.grid(which='both')
      plt.legend()
  
      print(&quot;Maximum of CDF is: {:2.1f}&quot;.format(the_cdf[-1]))
  
      plt.show()
  </code></pre>
  <p>The generated image is:
  <a href="https://i.stack.imgur.com/rZ8r0.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/rZ8r0.png" alt="enter image description here" /></a></p>
  <p>The output is:</p>
  <pre class="lang-sh prettyprint-override"><code>PDF calc time: 0.0010
  CDF calc time: 0.0010
  RVS calc time: 11.1120
  Maximum of CDF is: 1.0
  </code></pre>
  <p>The time computing the RVS method is too slow in my approach.</p>