• [x] Estimate statistical uncertainty of p-value from permutation tests
• [x] Allow users to pass a target precision to the permutation test instead of asking for the number of bootstrap samples
• [x] Generic permutation test for user-defined test statistics
• [x] allow sample arguments in resample, bootstrap, etc.
• [x] Include USP test (see #89)
• [x] Document changes and write migration guide
• [x] Fix the iterative algorithm to find the required number of permutation samples for usp and same_population

As pointed out in the Wikipedia, we should give a confidence interval for the computed pvalue in permutation tests.

With this one could also implement an automatic mode, where the number of required permutations is determined automatically that is needed to determine with near certainty whether the pvalue is below or above a given threshold α.

The pvalue confidence interval can be added as another field to the PermutationResult class. In order to use the automatic mode, the user has to provide alpha instead of size when the test is called, so the option for the test should probably be called size_or_alpha.

@HDembinski Couple edits to setup.py before releasing the new version (I can also make you co-owner on PyPI). Let me know if there are any other loose ends you can think of that need to be taken care of before releasing (it's been two years since the last time I did this and need to refresh my memory).

The speed of the functions jackknife and bootstrap could be greatly enhanced with numba. While it is possible to have two implements for each function, one requiring numba and one just requiring numpy, I think it would make sense to require numba. It used to be a pain to install, but nowadays one can just pip install it and of course it is also available from anaconda.

I know numba well enough to work on this. Another option (but less preferred), would be to write the implementations in C++ and wrap the code to Python with pybind11. I can do that as well, but then the package is not pure Python anymore and has to be compiled for every target platform. This is a huge added burden, and the C++ code may not run faster than numba (I did some tests a while ago, where numba beat my C++ implemenation).

I think there are several tools which fit very nicely in a separate module called empirical.

utils.ecdf -> empirical.cdf utils.eqf -> empirical.quantile jackknife.empirical_influence -> empirical.influence

@dsaxton Sound good?

Most professional repos use a two-branch model with master and develop (Boost, Scikit-HEP...), I think we should switch to that.

develop:

• default branch, commits and PRs go there
• force updates allowed but should be kept to an absolute minimum
• forks/feature branches are based on develop
• gets rebased to master if master was updated for a release
• RTD "latest" build points to develop

master:

• never force-updated (can and should be turned off in the settings)
• only updated on making a release
• on release, merge develop into master, update version number and changelog, make release
• special CI jobs that check the release are only triggered on updates to master
• RTD "stable" build points to master and is shown by default

The two main advantages of this scheme are that

• one can have special CI jobs that check releases, which are only triggered on updates to master
• one can have RTD point to the latest release automatically (using the "stable" build)

As mentioned in #34, we need a way to allow computation from pre-calculated replicates.

My original idea was to reuse the existing interface to do this, but documentation and argument types then become a bit ugly, which was also @dsaxton 's concern.

We need to resolve this before publishing release 1.0, in case it has implications for our interface overhaul (not necessarily the case, but it could have).

To address the points raised, I drafted a solution here where the interface of the bias function in the jackknife module is left as is. In addition, there is a bias_from_precalculated (name to be refined), which accepts the pre-calculated replicates. Internally, bias calls bias_from_precalculated, of course.

@dsaxton Would that be a way to go for all functions? We need to introduce X_from_precalculated then for

• bias
• bias_corrected
• variance
• confidence_interval

in both the jackknife and bootstrap modules.

This would be an acceptable solution for me, but whenever I see a common prefix/suffix, I am thinking of namespaces. I think it would be more organized to put these in a separate module, so that one can do

from resample.jackknife.precalculated import bias # the version in which you pass `theta` and `resampled_thetas`

or

from resample.jackknife import bias # the version in which you pass `fn` and `sample`

I think to make this work, we need to make jackknife and bootstrap into sub-packages, which then can have a sub-module precalculated. The directory structure would look like this.

resample
  __init__.py
  jackknife
     __init__.py
     precalculated.py
  bootstrap
     __init__.py
     precalculated.py

I took the liberty to submit an abstract to PyHEP 2020 virtual workshop on Python data analysis in high energy physics:

https://indico.cern.ch/event/882824/overview

I think this is a good place to advertise resample and what can be done with it. The workshop is very popular this year, so I went directly ahead in submitting the abstract without asking, but I can withdraw it if you don't think this is a good idea.

The presentation will mainly consist of a Jupyter notebook which show-cases the capabilities of resample, there will be some general information, too, like who is working on it, status, plans, etc.

tl;dr: I am a very experienced developer and a bootstrap expert, care to join forces?

I am the developer of Boost.Histogram and the maintainer of iminuit. I am a core member of the Scikit-HEP community. I am a senior developer in C++ and Python and a statistical expert. I read Efron's book "The introduction to the bootstrap" and a few others on resampling methods and think that resampling has a huge potential for my field (particle physics).

I think that we need a high-quality resampling package in Python, which ideally should become part of scipy at some point. I am privately using the balanced bootstrap since many years. I have an implementation of the balanced bootstrap in my tool library pyik. I would like to move away from my own implementation, however, because that is not sustainable. I would be happy to support one of the projects that focus entirely on resampling, such as yours. Your package seems to be the most advanced in terms of quality and completeness from what I have seen on PyPI. Are you interested in collaboration?

Best regards!

Closes https://github.com/dsaxton/resample/issues/37

I'm not 100% positive this will work, but it seems it should prevent it from using Python 2 for the build.

https://docs.readthedocs.io/en/stable/config-file/v2.html

I added a warning to the bias computation with the bootstrap. Wording should be further improved for clarity.

I discovered that the bootstrap cannot detect biases of the kind (unbiased estimate + constant / N), where N is the number of observations, which the jackknife does detect and remove exactly. @dsaxton You can see it for yourself if you apply the bootstrap bias computation on the test function def bad_mean(x): return (np.sum(x) + 2) / len(x) that I used to test the jackknife. The bootstrap bias considers this function unbiased.

Somehow the meaning of the two bias computations is different. The jackknife computes the bias of the estimator on the finite sample with respect to the asymptotic limit. The bootstrap bias seems to compute something different.

I am not sure if this question might be better placed on StackOverflow, but I guess it could be also a documentation issue so I guess it could also fit here. Besides, the code below might also help to explain what I meant with issue #112, so I hope, I can post it here on GitHub. I fitted a multinomial logistic regression to my data. Now I would like to get 95%-CI intervals for my beta-coefficients and additional p-values for each of them (and for each class, since with multinomial LR you get coefficients for each classes). How would I do that with your package?

Here's some example code (as written above, besides from this issue, it might also help to explain what I described in #112):

import numpy as np
from sklearn.datasets import load_iris

from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import Pipeline

from resample.bootstrap import bootstrap

# set random number generator
rng = np.random.RandomState(42)

# load data
X,y = load_iris(return_X_y=True)

# implement Pipeline
# center data and fit multinomial logistic regression
pipe = Pipeline([('scaler',StandardScaler()),
                  ('lr',LogisticRegression())])


pipe.fit(X,y)

# prepare data for bootstrapping with resample package
# NOTE: This is a little tedious because we have to concatenate X an y first
# to one array only to separate them again a few steps later within the function
A = np.concatenate((y.reshape(-1,1),X),axis=1)

# NOTE: Currently, the resample package does not offer to provide any further
# arguments that are passed to the function besides A. So we have to rely on 
# global scope (which might be a little bit unsafe?)
def fit_mlr(A):
    
    X = A[:,1:]
    y = A[:,0]
    pipe.fit(X,y)

    return {'coef': pipe._final_estimator.coef_,
            'intercept': pipe._final_estimator.coef_}

boot_coef = bootstrap(fit_mlr,sample=A,size=10)

I worked out how to iterate the permutation tests until a precision target is reached.

It would be great to implement such a functionality also for the functions

• bootstrap.bias
• bootstrap.bias_corrected
• bootstrap.variance
• bootstrap.confidence_interval

This means adding keywords precision and max_size to the functions and to deprecate size (which would act like max_size=size, precision=0). We cannot do this for bootstrap.bootstrap, because we don't know what the user is computing.

With the keyword return_error (default False) we can optionally return the calculated uncertainty in a backward compatible way.

One of the first questions I got after the presentation on resample at PyHEP was about parallelization.

In principle, resampling methods are perfectly parallelizable, assuming that fn is pure (has no side-effects). That is generally a reasonable assumption. In Python, there are many ways to parallelize, you may want to parallelize on your own cores, or on some cluster of computers, or on the cloud. Therefore, offering direct access to resample is good, because it allows the user to user to chose their parallelization scheme.

For the simple common cases, however, we may want to offer a threads option to our methods, which compute fn on the replicas using threads number of threads on the current computer, to better utilize common multi-core processors. This would an option for the functions bootstrap and jackknife and those that build on them, e.g. bias and variance etc. @dsaxton What do you think?

The jackknife-after-bootstrap method, as described in Efron and Tibshirani's book, is a clever way to compute an uncertainty for a bootstrap estimate, without computing additional replicas. It needs a bit of additional book-keeping, so it does not come for free, but it is a vast improvement over doing a full jackknife after the bootstrap.

We could add this an keyword option in resample.bootstrap.bootstrap, or have a separate resample.bootstrap.jackknife_after_bootstrap function. I am leaning slightly towards the latter.

I started to use the library in practice and I found a caveat of our current design. Let's say I want to compute a bias correction and the variance of an estimator. If I naively call resample.jackknife.variance and resample.jackknife.bias_corrected, it computes the jackknife estimates twice (which is expensive). The interface should allow me to reuse precomputed jackknife estimates (I am talking about the jackknife but the same is true for the bootstrap).

I am not sure yet how to best achieve this. Here is my idea so far.

Currently, we have in resample.jackknife the signature def variance(fn, sample). It expects two mandatory arguments and I think that should not change. However, we could make it so that if one passes None for fn, then sample is interpreted as the precomputed replicas. This is not ambiguous, because fn is never None under normal circumstances.

This approach works for all jackknife tools, but resample.bootstrap.confidence_level adds further complications. More precisely, when the "student" and "bca" methods are used, the baseline idea does not work. The "student" method also needs fn(sample) in addition to the replicas, and "bca" also needs fn(sample) and jackknife replicas on top.

I think the basic idea can still work, if we make the call to confidence_interval like this

thetas = bootstrap(my_fn, data)
theta = my_fn(data)
j_thetas = jackknife(my_fn, data)
confidence_interval(None, thetas, ci_method="percentile") # ok, works
confidence_interval(None, (thetas, theta), ci_method="student") # ok, additional information passed as tuple
confidence_invertal(None, (thetas, theta, j_thetas), ci_method="bca") # ok, same

Any thoughts?