code: def bovy_ars(domain,isDomainFinite,abcissae,hx,hpx,nsamples=1,
               hxparams=(),maxn=100):
      
      #First set-up the upper and lower hulls
      hull=setup_hull(domain,isDomainFinite,abcissae,hx,hpx,hxparams)
      #Then start  sampling: call sampleone repeatedly
      out= []
      nupdates= 0
      for ii in range(int(nsamples)):
          thissample, hull, nupdates= sampleone(hull,hx,hpx,domain,isDomainFinite,maxn,nupdates,hxparams)
          out.append(thissample)
      return out
docs: bovy_ars: Implementation of the Adaptive-Rejection Sampling
      algorithm by Gilks & Wild (1992): Adaptive Rejection Sampling
      for Gibbs Sampling, Applied Statistics, 41, 337
      Based on Wild & Gilks (1993), Algorithm AS 287: Adaptive Rejection
      Sampling from Log-concave Density Functions, Applied Statistics, 42, 701
  
      Input:
  
         domain          - [.,.] upper and lower limit to the domain
  
         isDomainFinite  - [.,.] is there a lower/upper limit to the domain?
  
         abcissae        - initial list of abcissae (must lie on either side of the peak in hx if the domain is unbounded
  
         hx              - function that evaluates h(x) = ln g(x)
  
         hpx             - function that evaluates hp(x) =  d h(x) / d x
  
         nsamples        - (optional) number of desired samples (default=1)
  
         hxparams        - (optional) a tuple of parameters for h(x) and h'(x)
  
         maxn            - (optional) maximum number of updates to the hull (default=100)
                           
      Output:
  
         list with nsamples of samples from exp(h(x))
  
      External dependencies:
  
         math
         scipy
         scipy.stats
  
      History:
         2009-05-21 - Written - Bovy (NYU)