biogeme.draws

Examples of use of several functions.

This is designed for programmers who need examples of use of the functions of the module. The examples are designed to illustrate the syntax. They do not correspond to any meaningful model.

author:

Michel Bierlaire

date:

Tue Nov 21 18:36:59 2023

import numpy as np
import pandas as pd
from biogeme.version import getText
import biogeme.draws as dr

Version of Biogeme.

print(getText())

biogeme 3.2.13 [2023-12-23]
Home page: http://biogeme.epfl.ch
Submit questions to https://groups.google.com/d/forum/biogeme
Michel Bierlaire, Transport and Mobility Laboratory, Ecole Polytechnique Fédérale de Lausanne (EPFL)

We set the seed so that the outcome of random operations is always the same.

np.random.seed(90267)

Uniform draws

Uniform [0,1]. The output is transformed into a data frame just for the display.

draws = dr.getUniform(sample_size=3, number_of_draws=10, symmetric=False)
pd.DataFrame(draws)

	0	1	2	3	4	5	6	7	8	9
0	0.406996	0.952091	0.393300	0.422424	0.392314	0.656022	0.071824	0.282152	0.294621	0.294378
1	0.719719	0.976471	0.860015	0.637804	0.097748	0.215859	0.901997	0.259457	0.246732	0.608131
2	0.914063	0.253449	0.299638	0.023243	0.008511	0.309165	0.172084	0.848053	0.775280	0.234141

draws = dr.getUniform(sample_size=3, number_of_draws=10, symmetric=True)
pd.DataFrame(draws)

	0	1	2	3	4	5	6	7	8	9
0	0.541595	-0.613439	0.410899	0.656565	-0.377678	0.978805	-0.199249	-0.499704	-0.941578	0.433283
1	-0.395955	0.442475	0.265263	0.511135	-0.507459	0.738372	-0.868453	-0.620988	0.553554	-0.213905
2	-0.797249	0.439322	0.127049	0.979762	0.609904	-0.430471	0.380190	0.536699	0.022698	-0.761527

LatinHypercube: the Modified Latin Hypercube Sampling (MLHS, Hess et al, 2006) provides U[0,1] draws from a perturbed grid, designed for Monte-Carlo integration.

latin_hypercube = dr.getLatinHypercubeDraws(sample_size=3, number_of_draws=10)
pd.DataFrame(latin_hypercube)

	0	1	2	3	4	5	6	7	8	9
0	0.325211	0.630001	0.177462	0.530946	0.417055	0.825831	0.782322	0.088749	0.155039	0.539467
1	0.986429	0.028337	0.585812	0.846593	0.661135	0.925723	0.397986	0.350411	0.955311	0.702116
2	0.744942	0.884753	0.046556	0.131465	0.463126	0.289113	0.685699	0.251115	0.204371	0.481789

The same method can be used to generate draws from U[-1,1]

latin_hypercube = dr.getLatinHypercubeDraws(
    sample_size=5, number_of_draws=10, symmetric=True
)
pd.DataFrame(latin_hypercube)

	0	1	2	3	4	5	6	7	8	9
0	-0.487827	0.946583	-0.290890	-0.997732	0.685466	-0.929016	0.815377	0.060902	-0.703703	-0.660290
1	-0.119771	0.903854	0.518566	-0.175317	-0.593672	-0.476686	-0.249879	0.014049	0.592416	0.333544
2	-0.536467	-0.602720	0.421731	0.258026	-0.369102	-0.215673	0.960123	-0.406895	0.652804	0.186110
3	0.149259	-0.053414	-0.356505	0.636484	0.521282	-0.866888	-0.784501	0.312416	-0.747227	0.112625
4	-0.914200	0.238546	0.738524	-0.027154	0.798928	-0.134140	-0.800939	0.460820	0.362292	0.852471

The user can provide her own series of U[0,1] draws.

my_unif = np.random.uniform(size=30)
pd.DataFrame(my_unif)

	0
0	0.061760
1	0.201607
2	0.779832
3	0.602527
4	0.696678
5	0.764274
6	0.568655
7	0.371194
8	0.233237
9	0.023673
10	0.339896
11	0.060995
12	0.221347
13	0.476015
14	0.866055
15	0.350926
16	0.644284
17	0.904649
18	0.857553
19	0.132332
20	0.408866
21	0.415191
22	0.502339
23	0.498379
24	0.612381
25	0.675059
26	0.349810
27	0.033002
28	0.145176
29	0.237861

latin_hypercube = dr.getLatinHypercubeDraws(
    sample_size=3, number_of_draws=10, symmetric=False, uniformNumbers=my_unif
)
pd.DataFrame(latin_hypercube)

	0	1	2	3	4	5	6	7	8	9
0	0.783279	0.092661	0.407378	0.274441	0.628585	0.901100	0.596822	0.300789	0.245706	0.218955
1	0.002059	0.750078	0.554809	0.156556	0.974595	0.192142	0.820413	0.040054	0.680296	0.938173
2	0.511698	0.449200	0.637744	0.120084	0.855835	0.878327	0.368700	0.713840	0.344663	0.495535

The uniform draws can also be arranged in a two-dimension array

my_unif = dr.getUniform(sample_size=3, number_of_draws=10)
pd.DataFrame(my_unif)

	0	1	2	3	4	5	6	7	8	9
0	0.429347	0.603875	0.896565	0.754603	0.769762	0.451352	0.832002	0.224726	0.904498	0.097305
1	0.091637	0.288086	0.133953	0.360001	0.793237	0.569249	0.338901	0.982264	0.660052	0.075265
2	0.336336	0.374114	0.464009	0.227873	0.008158	0.286625	0.934171	0.109902	0.945957	0.756240

latin_hypercube = dr.getLatinHypercubeDraws(
    sample_size=3, number_of_draws=10, uniformNumbers=my_unif
)
pd.DataFrame(latin_hypercube)

	0	1	2	3	4	5	6	7	8	9
0	0.181712	0.800272	0.903663	0.897806	0.493108	0.518975	0.774262	0.544630	0.599409	0.125153
1	0.677878	0.991875	0.748800	0.336388	0.445333	0.404465	0.712470	0.158992	0.014312	0.635842
2	0.096552	0.296817	0.053463	0.240824	0.303244	0.842887	0.227733	0.376270	0.964865	0.622002

Halton draws

One Halton sequence.

halton = dr.getHaltonDraws(sample_size=2, number_of_draws=10, base=3)
pd.DataFrame(halton)

	0	1	2	3	4	5	6	7	8	9
0	0.333333	0.666667	0.111111	0.444444	0.777778	0.222222	0.555556	0.888889	0.037037	0.370370
1	0.703704	0.148148	0.481481	0.814815	0.259259	0.592593	0.925926	0.074074	0.407407	0.740741

Several Halton sequences.

halton = dr.getHaltonDraws(sample_size=3, number_of_draws=10)
pd.DataFrame(halton)

	0	1	2	3	4	5	6	7	8	9
0	0.50000	0.25000	0.75000	0.12500	0.62500	0.37500	0.87500	0.06250	0.56250	0.31250
1	0.81250	0.18750	0.68750	0.43750	0.93750	0.03125	0.53125	0.28125	0.78125	0.15625
2	0.65625	0.40625	0.90625	0.09375	0.59375	0.34375	0.84375	0.21875	0.71875	0.46875

Shuffled Halton sequences.

halton = dr.getHaltonDraws(sample_size=3, number_of_draws=10, shuffled=True)
pd.DataFrame(halton)

	0	1	2	3	4	5	6	7	8	9
0	0.25000	0.12500	0.28125	0.46875	0.53125	0.59375	0.09375	0.31250	0.81250	0.21875
1	0.71875	0.78125	0.50000	0.37500	0.93750	0.62500	0.87500	0.90625	0.40625	0.65625
2	0.84375	0.18750	0.43750	0.56250	0.34375	0.03125	0.68750	0.15625	0.06250	0.75000

The above sequences were generated using the default base: 2. It is possible to generate sequences using different prime numbers.

halton = dr.getHaltonDraws(sample_size=1, number_of_draws=10, base=3)
pd.DataFrame(halton)

	0	1	2	3	4	5	6	7	8	9
0	0.333333	0.666667	0.111111	0.444444	0.777778	0.222222	0.555556	0.888889	0.037037	0.37037

It is also possible to skip the first items of the sequence. This is desirable in the context of Monte-Carlo integration.

halton = dr.getHaltonDraws(sample_size=1, number_of_draws=10, base=3, skip=10)
pd.DataFrame(halton)

	0	1	2	3	4	5	6	7	8	9
0	0.703704	0.148148	0.481481	0.814815	0.259259	0.592593	0.925926	0.074074	0.407407	0.740741

Antithetic draws

Antithetic draws can be generated from any function generating uniform draws.

draws = dr.getAntithetic(dr.getUniform, sample_size=3, number_of_draws=10)
pd.DataFrame(draws)

	0	1	2	3	4	5	6	7	8	9
0	0.422126	0.471794	0.148715	0.738601	0.041052	0.577874	0.528206	0.851285	0.261399	0.958948
1	0.512323	0.765719	0.501317	0.146468	0.558539	0.487677	0.234281	0.498683	0.853532	0.441461
2	0.417705	0.778921	0.752020	0.043203	0.255528	0.582295	0.221079	0.247980	0.956797	0.744472

Antithetic MLHS

draws = dr.getAntithetic(dr.getLatinHypercubeDraws, sample_size=3, number_of_draws=10)
pd.DataFrame(draws)

	0	1	2	3	4	5	6	7	8	9
0	0.078073	0.313795	0.158663	0.015185	0.421812	0.921927	0.686205	0.841337	0.984815	0.578188
1	0.478287	0.769973	0.244302	0.612203	0.906746	0.521713	0.230027	0.755698	0.387797	0.093254
2	0.949357	0.561688	0.722550	0.350440	0.807525	0.050643	0.438312	0.277450	0.649560	0.192475

Antithetic Halton.

draws = dr.getAntithetic(dr.getHaltonDraws, sample_size=1, number_of_draws=10)
pd.DataFrame(draws)

	0	1	2	3	4	5	6	7	8	9
0	0.5	0.25	0.75	0.125	0.625	0.5	0.75	0.25	0.875	0.375

As antithetic Halton draws may be correlated, it is a good idea to skip the first draws.

def uniform_halton(sample_size: int, number_of_draws: int) -> np.ndarray:
    """Function generating uniform draws for the antithetic draws"""
    return dr.getHaltonDraws(number_of_draws, sample_size, skip=100)

draws = dr.getAntithetic(uniform_halton, sample_size=3, number_of_draws=10)
pd.DataFrame(draws)

	0	1	2	3	4	5
0	0.648438	0.398438	0.898438	0.351562	0.601562	0.101562
1	0.085938	0.585938	0.335938	0.914062	0.414062	0.664062
2	0.835938	0.210938	0.710938	0.164062	0.789062	0.289062
3	0.460938	0.960938	0.054688	0.539062	0.039062	0.945312
4	0.554688	0.304688	0.804688	0.445312	0.695312	0.195312

Normal draws

Generate pseudo-random numbers from a normal distribution N(0,1) using the Algorithm AS241 Appl. Statist. (1988) Vol. 37, No. 3 by Wichura

draws = dr.getNormalWichuraDraws(sample_size=3, number_of_draws=10)
pd.DataFrame(draws)

	0	1	2	3	4	5	6	7	8	9
0	0.468384	0.309445	-0.727051	1.477318	2.324986	-0.385240	2.417537	-0.898517	1.043511	-0.527998
1	1.928242	0.392652	-1.698057	0.722431	0.231537	1.291789	-0.390652	0.757420	-0.111766	-0.253738
2	-1.628487	0.809681	-0.477319	-0.605592	-1.042455	2.172230	-1.201602	1.831441	-1.492897	0.730499

The antithetic version actually generates half of the draws and complete them with their antithetic version

draws = dr.getNormalWichuraDraws(sample_size=3, number_of_draws=10, antithetic=True)
pd.DataFrame(draws)

	0	1	2	3	4	5	6	7	8	9
0	-0.397450	-0.397369	1.096720	-0.516517	0.252910	0.397450	0.397369	-1.096720	0.516517	-0.252910
1	-0.318315	0.481563	0.992685	0.364958	0.601621	0.318315	-0.481563	-0.992685	-0.364958	-0.601621
2	-0.086190	0.780599	-0.178448	0.164590	-0.035292	0.086190	-0.780599	0.178448	-0.164590	0.035292

The user can provide her own series of U[0,1] draws. In this example, we use the MLHS procedure to generate these draws. Note that, if the antithetic version is used, only half of the requested draws must be provided.

my_unif = dr.getLatinHypercubeDraws(sample_size=3, number_of_draws=5)
pd.DataFrame(my_unif)

	0	1	2	3	4
0	0.730781	0.974142	0.544797	0.090617	0.654064
1	0.528956	0.221999	0.009099	0.420511	0.877654
2	0.154954	0.848118	0.760657	0.335545	0.326175

draws = dr.getNormalWichuraDraws(
    sample_size=3, number_of_draws=10, uniformNumbers=my_unif, antithetic=True
)
pd.DataFrame(draws)

	0	1	2	3	4	5	6	7	8	9
0	0.615176	1.945487	0.112528	-1.336963	0.396315	-0.615176	-1.945487	-0.112528	1.336963	-0.396315
1	0.072645	-0.765461	-2.361127	-0.200585	1.163339	-0.072645	0.765461	2.361127	0.200585	-1.163339
2	-1.015414	1.028396	0.708416	-0.424651	-0.450501	1.015414	-1.028396	-0.708416	0.424651	0.450501

The same with Halton draws.

my_unif = dr.getHaltonDraws(sample_size=2, number_of_draws=5, base=3, skip=10)
pd.DataFrame(my_unif)

	0	1	2	3	4
0	0.703704	0.148148	0.481481	0.814815	0.259259
1	0.592593	0.925926	0.074074	0.407407	0.740741

draws = dr.getNormalWichuraDraws(
    number_of_draws=10, sample_size=2, uniformNumbers=my_unif, antithetic=True
)
pd.DataFrame(draws)

	0	1	2	3	4	5	6	7	8	9
0	0.535083	-1.044409	-0.046436	0.895780	-0.645631	-0.535083	1.044409	0.046436	-0.895780	0.645631
1	0.234219	1.446104	-1.446104	-0.234219	0.645631	-0.234219	-1.446104	1.446104	0.234219	-0.645631