def monte_carlo_integration(distr_rvs, f, initial_chunk_size=100, chunk_size=10, eps=0.1):
    """
    Monte-Carlo numerical integration.
    :param distr_rvs: random sampler (takes size as `size` parameter and returns numpy array of samples
    :param f: function to be integrated
    :param initial_chunk_size: initial chunk size
    :param chunk_size: chunk size
    :param eps: precision
    :return: integral value
    """
    current_size = initial_chunk_size
    prev = 0.
    current = sum(map(f, distr_rvs(size=current_size))) / current_size
    i = 0
    while abs(current - prev) > eps:
        i += 1
        prev = current
        current = current * current_size / (current_size + chunk_size)
        current_size += chunk_size
        current += sum(map(f, distr_rvs(size=chunk_size))) / current_size
    return current