Various integration methods

Calculation of a simple integral using numerical integration and Monte-Carlo integration with various types of draws, including Halton draws base 13. It illustrates how to use draws that are not directly available in Biogeme.

Michel Bierlaire, EPFL Sat Jun 28 2025, 21:06:47

import numpy as np
import pandas as pd
from IPython.core.display_functions import display

from biogeme.biogeme import BIOGEME
from biogeme.database import Database
from biogeme.draws import RandomNumberGeneratorTuple, get_halton_draws
from biogeme.expressions import Draws, MonteCarlo, exp

We create a fake database with one entry, as it is required to store the draws.

df = pd.DataFrame()
df['FakeColumn'] = [1.0]
database = Database('fakeDatabase', df)

def halton13(sample_size: int, number_of_draws: int) -> np.array:
    """
    The user can define new draws. For example, Halton draws
    with base 13, skipping the first 10 draws.

    :param sample_size: number of observations for which draws must be
                       generated.
    :param number_of_draws: number of draws to generate.

    """
    return get_halton_draws(sample_size, number_of_draws, base=13, skip=10)

my_draws = {
    'HALTON13': RandomNumberGeneratorTuple(
        generator=halton13, description='Halton draws, base 13, skipping 10'
    )
}

Integrate with pseudo-random number generator.

integrand = exp(Draws('U', 'UNIFORM'))
simulated_integral = MonteCarlo(integrand)

Integrate with Halton draws, base 2.

integrand_halton = exp(Draws('U_halton', 'UNIFORM_HALTON2'))
simulated_integral_halton = MonteCarlo(integrand_halton)

Integrate with Halton draws, base 13.

integrand_halton13 = exp(Draws('U_halton13', 'HALTON13'))
simulated_integral_halton13 = MonteCarlo(integrand_halton13)

Integrate with MLHS.

integrand_mlhs = exp(Draws('U_mlhs', 'UNIFORM_MLHS'))
simulated_integral_mlhs = MonteCarlo(integrand_mlhs)

True value

true_integral = exp(1.0) - 1.0

Number of draws.

R = 2_000_000

sample_variance = (
    MonteCarlo(integrand * integrand) - simulated_integral * simulated_integral
)
stderr = (sample_variance / R) ** 0.5
error = simulated_integral - true_integral

sample_variance_halton = (
    MonteCarlo(integrand_halton * integrand_halton)
    - simulated_integral_halton * simulated_integral_halton
)
stderr_halton = (sample_variance_halton / R) ** 0.5
error_halton = simulated_integral_halton - true_integral

sampleVariance_halton13 = (
    MonteCarlo(integrand_halton13 * integrand_halton13)
    - simulated_integral_halton13 * simulated_integral_halton13
)
stderr_halton13 = (sampleVariance_halton13 / R) ** 0.5
error_halton13 = simulated_integral_halton13 - true_integral

sampleVariance_mlhs = (
    MonteCarlo(integrand_mlhs * integrand_mlhs)
    - simulated_integral_mlhs * simulated_integral_mlhs
)
stderr_mlhs = (sampleVariance_mlhs / R) ** 0.5
error_mlhs = simulated_integral_mlhs - true_integral

simulate = {
    'Analytical Integral': true_integral,
    'Simulated Integral': simulated_integral,
    'Sample variance   ': sample_variance,
    'Std Error         ': stderr,
    'Error             ': error,
    'Simulated Integral (Halton)': simulated_integral_halton,
    'Sample variance (Halton)   ': sample_variance_halton,
    'Std Error (Halton)         ': stderr_halton,
    'Error (Halton)             ': error_halton,
    'Simulated Integral (Halton13)': simulated_integral_halton13,
    'Sample variance (Halton13)   ': sampleVariance_halton13,
    'Std Error (Halton13)         ': stderr_halton13,
    'Error (Halton13)             ': error_halton13,
    'Simulated Integral (MLHS)': simulated_integral_mlhs,
    'Sample variance (MLHS)   ': sampleVariance_mlhs,
    'Std Error (MLHS)         ': stderr_mlhs,
    'Error (MLHS)             ': error_mlhs,
}

biosim = BIOGEME(
    database, simulate, random_number_generators=my_draws, number_of_draws=R
)
biosim.model_name = 'b02simple_integral'
results = biosim.simulate(the_beta_values={})
display(results)

   Analytical Integral  ...  Error (MLHS)
0             1.718282  ...              -2.157663e-10

[1 rows x 17 columns]

Reorganize the results.

print(f'Analytical integral: {results.iloc[0]["Analytical Integral"]:.6g}')
print(f"         \t{'Uniform':>15}\t{'Halton':>15}\t{'Halton13':>15}\t{'MLHS':>15}")
print(
    f'Simulated\t{results.iloc[0]["Simulated Integral"]:>15.6g}\t'
    f'{results.iloc[0]["Simulated Integral (Halton)"]:>15.6g}\t'
    f'{results.iloc[0]["Simulated Integral (Halton13)"]:>15.6g}\t'
    f'{results.iloc[0]["Simulated Integral (MLHS)"]:>15.6g}'
)
print(
    f'Sample var.\t{results.iloc[0]["Sample variance   "]:>15.6g}\t'
    f'{results.iloc[0]["Sample variance (Halton)   "]:>15.6g}\t'
    f'{results.iloc[0]["Sample variance (Halton13)   "]:>15.6g}\t'
    f'{results.iloc[0]["Sample variance (MLHS)   "]:>15.6g}'
)
print(
    f'Std error\t{results.iloc[0]["Std Error         "]:>15.6g}\t'
    f'{results.iloc[0]["Std Error (Halton)         "]:>15.6g}\t'
    f'{results.iloc[0]["Std Error (Halton13)         "]:>15.6g}\t'
    f'{results.iloc[0]["Std Error (MLHS)         "]:>15.6g}'
)
print(
    f'Error\t\t{results.iloc[0]["Error             "]:>15.6g}\t'
    f'{results.iloc[0]["Error (Halton)             "]:>15.6g}\t'
    f'{results.iloc[0]["Error (Halton13)             "]:>15.6g}\t'
    f'{results.iloc[0]["Error (MLHS)             "]:>15.6g}'
)

Analytical integral: 1.71828
                        Uniform          Halton        Halton13            MLHS
Simulated               1.71811         1.71828         1.71828         1.71828
Sample var.            0.241925        0.242035        0.242034        0.242036
Std error           0.000347797     0.000347876     0.000347875     0.000347876
Error              -0.000175252    -1.77867e-06    -4.40884e-06    -2.15766e-10