Vanilla Monte Carlo Integration

vanilla_monte_carlo.py

def MonteCarloIntegral(f, N, dim, fractional_error_tolerance=0.05):
        
    relative_error = 1.
    
    function_evaluations = list()
    while relative_error > fractional_error_tolerance:
        for _ in range(int(N)):
            F = f(np.random.rand(dim))
            function_evaluations.append(F)
            
        n_evaluations = len(function_evaluations)

        integral_estimate = np.mean(function_evaluations)
        variance_estimate = np.mean((np.array(function_evaluations) - integral_estimate)**2)
        error_estimate = np.sqrt(variance_estimate / n_evaluations)

        relative_error = error_estimate / integral_estimate

        
        # Calculate an estimate of how many extra points we would need to reach
        # our target relative error.
        #
        # var/N = (err)**2 and var/N_new = (new_err)**2 so:
        # N_new = (err/new_err)**2 * N
        # 
        # but we've already done N so the next iteration the requires: 
        # N * ((err/new_err)**2 - 1)
        # plus 1 to avoid case were we get stuck into an inf loop because N=0
        N = n_evaluations * ((relative_error / fractional_error_tolerance)**2 - 1) + 1
    
    print("Total evaluations:", n_evaluations)
    
    return (integral_estimate, variance_estimate, error_estimate)