Add DEMetropolisZ stepper

[michaelosthege]

The DE-MCMC-Z is a variant of DE-MCMC that uses its own trace/history to make new proposals instead of using a population of chains.

Compared to DEMetropolis this has several advantages:

• no need to step in parallel
• no between-process communication when running with cores > 1
• the "population" grows by every iteration -> works with chains=1 already

In the current implementation the history is kept in the instance of the stepper even though the same information is already present in the trace. But memory is cheap compared to inter-process communication & having to fiddle with sampling.py,

I'm working on a benchmarking/comparison notebook and will also do unit tests..

[codecov]

Codecov Report

Merging #3784 into master will increase coverage by 28.21%.
The diff coverage is 98.78%.

@@             Coverage Diff             @@
##           master    #3784       +/-   ##
===========================================
+ Coverage   62.57%   90.78%   +28.21%     
===========================================
  Files         133      133               
  Lines       20507    20537       +30     
===========================================
+ Hits        12832    18645     +5813     
+ Misses       7675     1892     -5783

Impacted Files	Coverage Δ
pymc3/tests/test_text_backend.py	100% <ø> (+100%)	⬆️
pymc3/tests/test_hmc.py	97.82% <ø> (+97.82%)	⬆️
pymc3/backends/text.py	96.07% <ø> (+74.5%)	⬆️
pymc3/distributions/dist_math.py	91.5% <ø> (+36.67%)	⬆️
pymc3/tests/test_quadpotential.py	95.45% <ø> (ø)	⬆️
pymc3/tests/helpers.py	56.6% <ø> (ø)	⬆️
pymc3/distributions/multivariate.py	78.93% <ø> (+9.57%)	⬆️
pymc3/tests/test_ndarray_backend.py	100% <ø> (+100%)	⬆️
pymc3/tests/test_model_func.py	100% <ø> (+100%)	⬆️
pymc3/tests/test_glm.py	100% <ø> (+100%)	⬆️
... and 250 more

[michaelosthege]

Here are some plots to better understand the strenghts and weaknesses of DEMetropolisZ:

Tuning

These tests are done with D=10, N_tune=10000, N_draws=10000, tune='lambda' and this model:

def get_mvnormal_model(D:int) -> pymc3.Model:
    true_mu = numpy.zeros(D)
    true_cov = numpy.eye(D)
    true_cov[:5,:5] = numpy.array([
        [1, 0.5, 0, 0, 0],
        [0.5, 2, 2, 0, 0],
        [0, 2, 3, 0, 0],
        [0, 0, 0, 4, 4],
        [0, 0, 0, 4, 5],
    ])

    with pymc3.Model() as pmodel:
        x = pymc3.MvNormal('x', mu=true_mu, cov=true_cov, shape=(D,))
    
    true_samples = x.random(size=1000)
    truth_id = arviz.data.convert_to_inference_data(true_samples[numpy.newaxis,:], group='random')
    return pmodel, truth_id

The tune_drop_fraction parameter helps to increase effective sample sizes by discarding parts of the tuning history.

When the entire tuning history is removed, the autocorrelation goes up:

Similarly, the acceptance rate jumps up:

Tuning lambda helps the initial convergence (same run as above). The chains were started at [5.0]*N_chains to put them outside of the typical set.

Efficiency compared to DEMetropolis

Compared to DEMetropolis, the DEMetropolisZ achieves the same effective sample sizes independent of the number of chains. At N_chains=2*D, DEMetropolis almost catches up:

I ran all experiments with cores=1 and there was no significant runtime difference.

[junpenglao]

Thanks @michaelosthege - I think comments like #3784 (comment) would be best in our wiki, so we can link it in the docstring if necessary.

[michaelosthege]

Thanks @michaelosthege - I think comments like #3784 (comment) would be best in our wiki, so we can link it in the docstring if necessary.

I have two or three notebooks for running these experiments that I would like to contribute in some way.
Maybe the "Deep Dive" section in the Tutorials ?

Also I should probably move the "pymc3.ode: Shapes and benchmarking" Example over to that section too.

[canyon289]

Very well written! Thanks for adding to PyMC3! Left some code style comments. Great work!!

[aloctavodia]

minor details