title: Speed up Metropolis--Hastings in Python
body: <p>I have some code that samples a posterior distribution using MCMC, specifically <a href="https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm" rel="nofollow noreferrer">Metropolis Hastings</a>. I use scipy to generate random samples:</p>
  
  <pre><code>import numpy as np
  from scipy import stats
  
  def get_samples(n):
      """
      Generate and return a randomly sampled posterior.
  
      For simplicity, Prior is fixed as Beta(a=2,b=5), Likelihood is fixed as Normal(0,2)
  
      :type n: int
      :param n: number of iterations
  
      :rtype: numpy.ndarray
      """
      x_t = stats.uniform(0,1).rvs() # initial value
      posterior = np.zeros((n,))
      for t in range(n):
          x_prime = stats.norm(loc=x_t).rvs() # candidate
          p1 = stats.beta(a=2,b=5).pdf(x_prime)*stats.norm(loc=0,scale=2).pdf(x_prime) # prior * likelihood 
          p2 = stats.beta(a=2,b=5).pdf(x_t)*stats.norm(loc=0,scale=2).pdf(x_t) # prior * likelihood 
          alpha = p1/p2 # ratio
          u = stats.uniform(0,1).rvs() # random uniform
          if u &lt;= alpha:
              x_t = x_prime # accept
              posterior[t] = x_t
          elif u &gt; alpha:
              x_t = x_t # reject
      posterior = posterior[np.where(posterior &gt; 0)] # get rid of initial zeros that don't contribute to distribution
      return posterior
  </code></pre>
  
  <p>Generally, I try to avoid using explicit for loops in python - I would try to generate everything using pure numpy. For the case of this algorithm however, a for loop with if statements is unavoidable. Therefore, the code is quite slow. When I profile my code, it is spending most time within the for loop (obviously), and more specifically, the slowest parts are in generating the random numbers; <code>stats.beta().pdf()</code> and <code>stats.norm().pdf()</code>.</p>
  
  <p>Sometimes I use <a href="http://numba.pydata.org/" rel="nofollow noreferrer">numba</a> to speed up my code for matrix operations. While numba is compatible with some numpy operations, generating random numbers is not one of them. Numba has a <a href="https://numba.pydata.org/numba-doc/dev/cuda/random.html" rel="nofollow noreferrer">cuda rng</a>, but this is limited to normal and uniform distributions.</p>
  
  <p>My question is, is there a way to significantly speed up the code above using some sort of random sampling of various distributions compatible with numba?</p>
  
  <p>We don't have to limit ourselves to numba, but it is the only easy to use optimizer that I know of. More generally, I am looking for ways to speed up random sampling for various distributions (beta, gamma, poisson) <strong>within a for loop</strong> in python.</p>