ENH: add stats.qmc module with quasi Monte Carlo functionality

[tupui]

Add scipy.stats.qmc.

Reference issue

Closes #9695

What does this implement/fix?

Provide a set of functions to create and assess quasi-Monte Carlo Design of Experiments.

It provides a generic class scipy.stats.qmc.QMCEngine which defines a QMC engine/sampler. An engine is state aware: it can be continued, advanced and reseted. 4 base samplers are available:

• scipy.stats.qmc.Sobol the well known Sobol' low discrepancy sequence. Several warnings have been added to guide the user into properly using this sampler. The sequence is scrambled by default.
• scipy.stats.qmc.Halton: Halton low discrepancy sequence. The sequence is scrambled by default.
• scipy.stats.qmc.LatinHypercube: plain LHS design.
• scipy.stats.qmc.OrthogonalLatinHypercube: Orthogonal version of LHS which is more uniform than plain LHS.

And 2 special samplers are available:

• scipy.stats.qmc.MultinomialQMC: sampling from a multinomial distribution using any of the base scipy.stats.qmc.QMCEngine.
• scipy.stats.qmc.MultivariateNormalQMC: sampling from a multivariate Normal using any of the base scipy.stats.qmc.QMCEngine.

The module also provide the following helpers:

• scipy.stats.qmc.discrepancy: assess the quality of a set of points in terms of space coverage.
• scipy.stats.qmc.update_discrepancy: can be used in an optimization loop to construct a good set of points.
• scipy.stats.qmc.scale: easily scale a set of points from (to) the unit interval to (from) a given range.

The module quality has been assessed by specialists from the domain (see the long discussion bellow for more details). Following are some convergence plots of Sobol' and Halton.

• gist for Halton: https://gist.github.com/tupui/ee6e685219c06c429338569e7ee7ff2a
• gist for Sobol': https://gist.github.com/tupui/fb6e219b1dd2316b7498ebce231bfff5

[rgommers]

Note: the scope of this PR is under discussion in gh-9695.

[tupui]

@rgommers First proposal is here. I still have a doc issue but I am using the same docstring as in scipy.optimize._differentialevolution.
I did not yet updated scipy.optimize._differentialevolution to use the new scipy.stats.qmc.orthogonal_latin_hypercube but appart for possible test fail there is no blocking point here.

[rgommers]

Just some initial comments. Scope discussion would be good to finish first before reviewing in more detail.

[tupui]

@tylerjereddy Hey I integrated your suggestions, thanks! Could you have another look?

[tupui]

@Balandat I finished to implement the changes we discussed.

• I moved Sobol to the file qmc (as well as MultivariateNormalQMC, NormalQMC, multinomial_qmc, discrepancy_star_L2). I had to change some cdef functions to cpdef for that. I think it's clearer like this.
• All the tests are updated.
• MultivariateNormalQMC, NormalQMC and multinomial_qmc accept an engine parameter. It falls back to Sobol.

I like the interface and I think this can be ready soon 😃. Let me know what you think.

[Balandat]

This looks great. FWIW I have a local commit that moves Sobol to use the new 21201-dimensional direction numbers from https://web.maths.unsw.edu.au/~fkuo/sobol/. The c-code generated from putting the numbers in as literals is obscenely large/inefficient, so what I did was to package and save the arrays as an .npz data file and load these in upon importing the Sobol generator. This works locally with fixing the path, the outstanding thing here is to make sure this is portable and the data file gets packaged with the extension and the path is resolved appropriately. Then the tests need to be modified to check fo the new direction numbers.

[Balandat]

Overall this looks good, couple of comments inline. There is a doc failure that needs to be fixed.

[tupui]

Overall this looks good, couple of comments inline. There is a doc failure that needs to be fixed.

I addressed the comments and fixed the doc (to be confirmed with tests). I plan to complete the docstrings with examples as with Halton.

Speaking of the doc. Is the current doc enough or should we add a bigger section somewhere about DoE, an introduction like for optimization? Do we want figures of DoE for instance? A section in https://docs.scipy.org/doc/scipy-1.4.1/reference/tutorial/stats.html ?

I have quite some material in case we need something a bit more theoretical (I have some chapters about DoE in my PhD thesis and articles I wrote, so license-wise fine).

cc @josef-pkt (I think you're the PoC for this part of the doc)

[tupui]

There is still a docstring warning about np.random.RandomState. It's in all seed param definition. But I am using here what's in the code of scipy.optimize (for instance).

I found another block used in rv_continuous. Maybe more up to date:

"""
seed : {None, int, `~np.random.RandomState`, `~np.random.Generator`}, optional
        This parameter defines the object to use for drawing random variates.
        If `seed` is `None` the `~np.random.RandomState` singleton is used.
        If `seed` is an int, a new ``RandomState`` instance is used, seeded
        with seed.
        If `seed` is already a ``RandomState`` or ``Generator`` instance,
        then that object is used.
        Default is None.
"""

[EDIT] fixed with this docstring

[tupui]

@Balandat I added quite some doc trying to mimic the style of scipy. IMO, we are good for a first version and this is ready for a review for integration 😃 .

[rgommers]

Did one more review pass through, resolved some comments that were addressed, and added stats.qmc to the list of API modules documented as public.

[rgommers]

Merged now, thank you @tupui and @Balandat, this is a really nice addition to SciPy! And thanks also to @mdhaber and @tylerjereddy for reviewing, and @ArtOwen, @fjhickernell, @rkern, @bashtage and others for contributing their expertise. This will be a highlight of the 1.7.0 release.

[rgommers]

Thanks for the patience on this one all!

[rgommers]

@tupui could you help with the following:

• Adding a description to the release notes (https://github.com/scipy/scipy/wiki/Release-note-entries-for-SciPy-1.7.0). I added it as a highlight, and added a subsection under scipy.stats improvements with TODO. A one paragraph summary would be great.
• Update the PR description with a short summary as well. The latest version of the convergence plot could be copied into that PR description too, it's quite nice and that makes it discoverable. Also, maybe that is useful for tutorial content?

[tupui]

Merged now, thank you @tupui and @Balandat, this is a really nice addition to SciPy! And thanks also to @mdhaber and @tylerjereddy for reviewing, and @ArtOwen, @fjhickernell, @rkern, @bashtage and others for contributing their expertise. This will be a highlight of the 1.7.0 release.

Thanks everyone for all the good work and your guidance during this process! I learned a lot during the process and I am grateful for this. I would also like to acknowledge Sergei Kucherenko who's been of great help via email.

For information, @Balandat and I will present this at the 13th International Conference on Monte Carlo Methods and Applications (MCM 2021). This for sure will help to spread the word about this new module.

[tupui]

@tupui could you help with the following:

• Adding a description to the release notes (https://github.com/scipy/scipy/wiki/Release-note-entries-for-SciPy-1.7.0). I added it as a highlight, and added a subsection under scipy.stats improvements with TODO. A one paragraph summary would be great.

• Update the PR description with a short summary as well. The latest version of the convergence plot could be copied into that PR description too, it's quite nice and that makes it discoverable. Also, maybe that is useful for tutorial content?

Sure I will open à PR with all that! I like the idea of the tutorial. I guess this would be another section of the existing stats tutorial.

[Balandat]

Yay! Don't think I've had >600 conversation items on a PR before, thanks eveyone for the thorough discussion!

[rgommers]

Sure I will open à PR with all that! I like the idea of the tutorial. I guess this would be another section of the existing stats tutorial.

Awesome. Yes indeed. At some point we may want to break this into separate pages for subsections, given that it is already very long. But for now I would just add to it.

Yay! Don't think I've had >600 conversation items on a PR before, thanks eveyone for the thorough discussion!

I know it was a painful delivery. We only get to add whole new submodules very rarely though:)

[Balandat]

I know it was a painful delivery. We only get to add whole new submodules very rarely though:)

Challenge accepted, next up: A full module.

[fjhickernell]

Thank you! Best, Fred Fred J. Hickernell Vice Provost for Research, research.iit.edu <https://research.iit.edu/> www.iit.edu/research/about/magazine <https://www.iit.edu/research/about/magazine> Professor, Department of Applied Mathematics Center for Interdisciplinary Scientific Computation (CISC) cos.iit.edu/cisc/ <https://cos.iit.edu/cisc/> Illinois Institute of Technology IIT Tower, Room 7C9-2, 10 W 35th St, Chicago, IL 60616 hickernell@iit.edu <mailto:hickernell@iit.edu>, +1 312 567 3470 <tel:+1%20312%20567%203470> Working remotely. Please contact me by email or phone.

…

On Jan 30, 2021, at 10:51 AM, Max Balandat ***@***.***> wrote: I know it was a painful delivery. We only get to add whole new submodules very rarely though:) Challenge accepted, next up: A full module. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#10844 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAGHS6VBPX4DIUXZLPRQJBTS4Q2HLANCNFSM4IYKCGMA>.

[tupui]

@rgommers I updated the description of the issue. Added some comments in the release notes (mostly the same text). And I created the 2 follow up PR I had. I will work on the tutorial now.