<p>There are several things you can do here.</p>
<ul>
<li>For fixed parameters <code>v_esc</code> and <code>v_0</code>, <code>n_0</code> is a constant, so it doesn't need to be calculated in the <code>pdf</code> method.</li>
<li>If you define only a PDF for a SciPy <code>rv_continuous</code> subclass, then the class's <code>rvs</code>, <code>mean</code>, and so on will be very slow, presumably because the method needs to integrate the PDF every time it needs to generate a random variate or calculate a statistic. If speed is at a premium, you will thus need to add to <code>maxwell_boltzmann_pdf</code> an <code>_rvs</code> method that uses its own sampler. (See also <a href="https://stackoverflow.com/questions/62516272">this question</a>.) One possible method is the <em>rejection sampling</em> method: Draw a number in a box until the box falls within the PDF. It works for any bounded PDF with a finite domain, as long as you know what the domain and bound are (the bound is the maximum value of <code>f</code> in the domain). See <a href="https://stackoverflow.com/questions/55724501">this question</a> for example Python code.</li>
<li>If you know the distribution's CDF, then there are some additional tricks. One of them is the relatively new <a href="https://arxiv.org/abs/2004.02339" rel="nofollow noreferrer"><em>k-vector sampling</em></a> method for sampling a continuous distribution. There are two phases: a setup phase and a sampling phase. The setup phase involves approximating the CDF's inverse via root finding, and the sampling phase uses this approximation to generate random numbers that follow the distribution in a very fast way without having to further evaluate the CDF. For a fixed distribution like this one, if you show me the CDF, I can precalculate the necessary data and the code needed to sample the distribution using that data. Essentially, the only non-trivial part of <em>k-vector sampling</em> is the root-finding step.</li>
<li>More information on sampling from an arbitrary distribution is found on my <a href="https://peteroupc.github.io/randomfunc.html#Random_Numbers_from_an_Arbitrary_Distribution" rel="nofollow noreferrer">sampling methods page</a>.</li>
</ul>
Getting random numbers from a truncated Maxwell-Boltzmann distribution
<p>I would like to generate random numbers using a truncated Maxwell-Boltzmann distribution. I know that scipy has built-in Maxwell random variables, but there is no truncated version of it (I am also aware of a truncated normal distribution, which is irrelevant here). I have tried to write my own random variables using rvs_continuous:</p>

<pre><code>import scipy.stats as st

class maxwell_boltzmann_pdf(st.rv_continuous):

    def _pdf(self,x):
        n_0 = np.power(np.pi,3/2)*np.square(v_0)*(v_0*erf(v_esc/v_0)-(2/np.sqrt(np.pi))*v_esc*np.exp(-np.square(v_esc/v_0)))
        return (1/n_0)*(4*np.pi*np.square(x))*np.exp(-np.square(x/v_0))*np.heaviside(v_esc-x,0)

maxwell_boltzmann_cv = maxwell_boltzmann_pdf(a=0, b=v_esc, name='maxwell_boltzmann_pdf')
</code></pre>

<p>This does exactly what I want, but it is way too slow for my purpose (I am doing Monte Carlo simulations), even if I draw all the random velocities outside of all the loops. I have also thought of using Inverse transform sampling method, but the inverse of the CDF does not have an analytic form and I will need to do a bisection for every number I draw. It would be great if there is a convenient way for me to generate random numbers from a truncated Maxwell-Boltzmann distribution with decent speed.</p>