title: Efficiently sample from arbitrary multivariate function
body: <p>I would like to sample from an arbitrary function in Python.</p>
  
  <p>In <a href="https://stackoverflow.com/questions/21100716/fast-arbitrary-distribution-random-sampling">Fast arbitrary distribution random sampling</a> it was stated that one could use inverse transform sampling and in <a href="https://stackoverflow.com/questions/4113307/pythonic-way-to-select-list-elements-with-different-probability">Pythonic way to select list elements with different probability</a> it was mentioned that one should use inverse cumulative distribution function. As far as I undestand those methods only work the univariate case. My function is multivariate though and too complex that any of the suggestions in <a href="https://stackoverflow.com/a/48676209/4533188">https://stackoverflow.com/a/48676209/4533188</a> would apply.</p>
  
  <p>Prinliminaries: My function is based on Rosenbrock's banana function, which value we can get the value of the function with</p>
  
  <pre><code>import scipy.optimize
  scipy.optimize.rosen([1.1,1.2])
  </code></pre>
  
  <p>(here <code>[1.1,1.2]</code> is the input vector) from scipy, see <a href="https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.optimize.rosen.html" rel="nofollow noreferrer">https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.optimize.rosen.html</a>.</p>
  
  <p>Here is what I came up with: I make a grid over my area of interest and calculate for each point the function value. Then I sort the resulting data frame by the value and make a cumulative sum. This way we get "slots" which have different sizes - points which have large function values have larger slots than points with small function values. Now we generate random values and look into which slot the random value falls into. The row of the data frame is our final sample.</p>
  
  <p>Here is the code:</p>
  
  <pre><code>import scipy.optimize
  from itertools import product
  from dfply import *
  
  nb_of_samples = 50
  nb_of_grid_points = 30
  
  rosen_data = pd.DataFrame(array([item for item in product(*[linspace(fm[0], fm[1], nb_of_grid_points) for fm in zip([-2,-2], [2,2])])]), columns=['x','y'])
  rosen_data['z'] = [np.exp(-scipy.optimize.rosen(row)**2/500) for index, row in rosen_data.iterrows()]
  rosen_data = rosen_data &gt;&gt; \
      arrange(X.z) &gt;&gt; \
      mutate(z_upperbound=cumsum(X.z)) &gt;&gt; \
      mutate(z_upperbound=X.z_upperbound/np.max(X.z_upperbound))
  value = np.random.sample(1)[0]
  
  def get_rosen_sample(value):
      return (rosen_data &gt;&gt; mask(X.z_upperbound &gt;= value) &gt;&gt; select(X.x, X.y)).iloc[0,]
  
  values = pd.DataFrame([get_rosen_sample(s) for s in np.random.sample(nb_of_samples)])
  </code></pre>
  
  <p>This works well, but I don't think it is very efficient. What would be a more efficient solution to my problem?</p>
  
  <p>I read that Markov chain Monte Carlo might helping, but here I am in over my head for now on how to do this in Python.</p>