Resample the output a function until the value convergences to given tolerance

README.md

Resample function until convergence

Python function that will rerun a function until its output convergences. The output must be numerical values. It can been an int, float, list, or array. Requires numpy.

This method runs the function in chunks (N=sample_size). It calculates the mean and standard deviation of the total set of samples after each chunk. If the mean and standard deviation of the total set is relatively close (see rtol) the last iteration, the sampling is assumed to have converged.

resample_to_convergence.py

import numpy as np

def resample_to_convergence(function, sample_size=50, rtol=10e-2, *args, **kwargs):
    """Repeat a calculation until the mean and standard deviations of the samples
    converge to a specified tolerance level (note that the tolerances are relative.)

    Parameters
    ----------
    function : callable
        Function to repeat until convergence.
    sample_size : int
        number of samples to draw before recalculating the statistics for convergence.
    rtol : float
        relative tolerance for convergence over sampling

    *args : position arguments for the callable function that is resampled
    **kwargs : keyword arguments for the callable function that is sampled.

    Returns
    -------
    samples : numpy.array
        array containing the samples. Each column is a sample. Rows are the different
        values from the output of the callable
    mean : numpy.array
        means of the samples
    std : numpy.array
        standard deviations of samples
    count : int
        number of samples needed to reach convergence test.

    """
    tests = (False,False)
    count = 0

    while False in tests:
        ############### Create a sample ###############
        init_sample = function(*args, **kwargs)
        # If init_sample is array
        try:
            # Allocate a 2d array with proper size.
            new_samples = np.empty((sample_size, len(init_sample)), dtype=float)

            # Populate the array with various samples
            new_samples[0, :] = init_sample
            for i in range(1,sample_size):
                new_samples[i, :] = function(*args, **kwargs)

        # If it is a value
        except TypeError:
            # Allocate memory in 1d array
            new_samples = np.empty(sample_size, float)
            new_samples[0] = init_sample

            for i in range(1,sample_size):
                new_samples[i] = function(*args, **kwargs)

        ############ Check for convergence ##############
        try:
            samples = np.concatenate((samples, new_samples), axis=0)
            # Calculate the statistics on the mean and sample.
            mean = np.mean(samples, axis=0)
            std = np.std(samples, ddof=1, axis=0)

            # Check values for
            check_mean = np.isclose(mean, old_mean, rtol=rtol)
            check_std = np.isclose(std, old_std, rtol=rtol)
            test1 = np.all(check_mean)
            test2 = np.all(check_std)
            tests = (test1, test2)

        # If this is the first iteration, don't check for convergence.
        except NameError:
            samples = new_samples
            mean = np.mean(samples, axis=0)
            std = np.std(samples, ddof=1, axis=0)

        old_mean = mean
        old_std = std

        count += sample_size

    return samples, mean, std, count