The cdf(method='exact') will still return the integrate.quad solution. I wasn't too happy with the cumsum implementation, but my idea was to prioritize speed first, and then have optional accuracy.

I was under the impression the density grid was uniform (my initial scan shows it being  build by a np.linspace call).