Draws: examples of use of each function

This webpage is for programmers who need examples of use of the functions of the package. The examples are designed to illustrate the syntax. They do not correspond to any meaningful model. For examples of models, visit biogeme.epfl.ch.

In [1]:
import datetime
print(datetime.datetime.now())

2019-12-29 21:16:55.723146
In [2]:
import biogeme.version as ver
print(ver.getText())

biogeme 3.2.5 [2019-12-29]
Version entirely written in Python
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)
In [3]:
import biogeme.draws as dr
import numpy as np
np.random.seed(90267)

Uniform [0,1]

In [4]:
draws = dr.getUniform(sampleSize=3,numberOfDraws=10,symmetric=False)
draws
Out[4]:
array([[0.4069955 , 0.95209084, 0.39330002, 0.42242449, 0.39231397,
        0.65602235, 0.07182417, 0.28215248, 0.29462105, 0.29437835],
       [0.71971885, 0.97647083, 0.86001524, 0.63780397, 0.09774763,
        0.21585889, 0.90199693, 0.25945665, 0.24673199, 0.60813096],
       [0.9140627 , 0.25344881, 0.29963774, 0.02324269, 0.00851069,
        0.3091653 , 0.17208375, 0.84805301, 0.77527991, 0.23414075]])
In [5]:
draws = dr.getUniform(sampleSize=3,numberOfDraws=10,symmetric=True)
draws
Out[5]:
array([[ 0.54159507, -0.61343854,  0.41089924,  0.65656456, -0.37767775,
         0.9788048 , -0.19924935, -0.49970433, -0.94157778,  0.43328274],
       [-0.39595535,  0.44247489,  0.26526338,  0.51113532, -0.50745861,
         0.73837224, -0.86845256, -0.62098784,  0.55355381, -0.21390505],
       [-0.7972486 ,  0.43932172,  0.12704909,  0.97976216,  0.60990425,
        -0.43047078,  0.3801898 ,  0.53669883,  0.02269837, -0.76152749]])

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.

In [6]:
latinHypercube = dr.getLatinHypercubeDraws(sampleSize=3,numberOfDraws=10)
latinHypercube
Out[6]:
array([[0.32521073, 0.63000128, 0.17746161, 0.53094562, 0.41705476,
        0.82583115, 0.78232205, 0.08874858, 0.15503933, 0.53946678],
       [0.9864292 , 0.02833746, 0.58581152, 0.84659261, 0.66113494,
        0.92572347, 0.39798632, 0.35041124, 0.9553105 , 0.70211636],
       [0.74494216, 0.88475329, 0.04655552, 0.13146474, 0.46312566,
        0.28911263, 0.68569908, 0.25111488, 0.204371  , 0.48178854]])

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

In [7]:
latinHypercube = dr.getLatinHypercubeDraws(sampleSize=5,numberOfDraws=10,symmetric=True)
latinHypercube
Out[7]:
array([[-0.48782677,  0.94658307, -0.29089021, -0.99773248,  0.68546628,
        -0.92901557,  0.81537697,  0.06090223, -0.70370281, -0.66028972],
       [-0.11977072,  0.90385423,  0.5185657 , -0.17531685, -0.59367191,
        -0.4766856 , -0.24987929,  0.01404931,  0.59241613,  0.33354383],
       [-0.53646676, -0.60271978,  0.42173134,  0.25802578, -0.36910209,
        -0.21567312,  0.9601231 , -0.40689483,  0.65280404,  0.18610977],
       [ 0.14925908, -0.0534143 , -0.35650472,  0.63648376,  0.52128198,
        -0.86688846, -0.78450064,  0.31241634, -0.74722651,  0.11262522],
       [-0.91420031,  0.23854628,  0.73852358, -0.02715404,  0.79892765,
        -0.13414034, -0.80093917,  0.46081992,  0.36229167,  0.85247066]])

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

In [8]:
myUnif = np.random.uniform(size=30)
myUnif
Out[8]:
array([0.06175992, 0.20160719, 0.77983176, 0.6025268 , 0.696678  ,
       0.76427402, 0.56865543, 0.37119385, 0.23323665, 0.0236732 ,
       0.33989596, 0.06099459, 0.2213469 , 0.47601499, 0.86605521,
       0.35092633, 0.64428447, 0.90464931, 0.85755345, 0.13233163,
       0.4088656 , 0.41519067, 0.50233919, 0.49837905, 0.61238122,
       0.67505865, 0.34980957, 0.03300237, 0.14517637, 0.23786104])
In [9]:
latinHypercube = dr.getLatinHypercubeDraws(sampleSize=3,numberOfDraws=10,symmetric=False,uniformNumbers=myUnif)
latinHypercube
Out[9]:
array([[0.7832793 , 0.09266106, 0.40737823, 0.27444122, 0.62858511,
        0.90110008, 0.59682164, 0.30078911, 0.24570646, 0.21895518],
       [0.00205866, 0.75007797, 0.55480948, 0.15655593, 0.97459537,
        0.19214247, 0.82041271, 0.04005357, 0.68029552, 0.93817255],
       [0.51169754, 0.4492005 , 0.63774439, 0.12008423, 0.85583529,
        0.87832699, 0.36869982, 0.71383969, 0.3446632 , 0.49553517]])

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

In [10]:
myUnif = dr.getUniform(sampleSize=3,numberOfDraws=10)
myUnif
Out[10]:
array([[0.42934713, 0.60387549, 0.8965652 , 0.75460304, 0.76976185,
        0.45135228, 0.83200157, 0.22472631, 0.90449758, 0.09730533],
       [0.09163669, 0.28808628, 0.13395278, 0.36000124, 0.79323743,
        0.56924873, 0.33890064, 0.98226408, 0.66005216, 0.07526522],
       [0.33633631, 0.37411379, 0.46400897, 0.22787268, 0.0081577 ,
        0.28662478, 0.93417098, 0.10990211, 0.94595662, 0.75624039]])
In [11]:
latinHypercube = dr.getLatinHypercubeDraws(sampleSize=3,numberOfDraws=10,uniformNumbers=myUnif)
latinHypercube
Out[11]:
array([[0.18171174, 0.80027192, 0.9036634 , 0.8978057 , 0.49310791,
        0.51897496, 0.77426242, 0.54463002, 0.5994088 , 0.12515343],
       [0.67787788, 0.99187468, 0.7488003 , 0.33638789, 0.44533337,
        0.40446509, 0.71247046, 0.15899206, 0.01431157, 0.63584217],
       [0.09655217, 0.29681659, 0.05346252, 0.24082421, 0.30324351,
        0.84288749, 0.22773339, 0.37626954, 0.96486522, 0.62200174]])

Halton draws

One Halton sequence

In [12]:
halton = dr.getHaltonDraws(sampleSize=2,numberOfDraws=10,base=3)
halton
Out[12]:
array([[0.33333333, 0.66666667, 0.11111111, 0.44444444, 0.77777778,
        0.22222222, 0.55555556, 0.88888889, 0.03703704, 0.37037037],
       [0.7037037 , 0.14814815, 0.48148148, 0.81481481, 0.25925926,
        0.59259259, 0.92592593, 0.07407407, 0.40740741, 0.74074074]])

Several Halton sequences

In [13]:
halton = dr.getHaltonDraws(sampleSize=3,numberOfDraws=10)
halton
Out[13]:
array([[0.5    , 0.25   , 0.75   , 0.125  , 0.625  , 0.375  , 0.875  ,
        0.0625 , 0.5625 , 0.3125 ],
       [0.8125 , 0.1875 , 0.6875 , 0.4375 , 0.9375 , 0.03125, 0.53125,
        0.28125, 0.78125, 0.15625],
       [0.65625, 0.40625, 0.90625, 0.09375, 0.59375, 0.34375, 0.84375,
        0.21875, 0.71875, 0.46875]])

Shuffled Halton sequences

In [14]:
halton = dr.getHaltonDraws(sampleSize=3,numberOfDraws=10,shuffled=True)
halton
Out[14]:
array([[0.25   , 0.125  , 0.28125, 0.46875, 0.53125, 0.59375, 0.09375,
        0.3125 , 0.8125 , 0.21875],
       [0.71875, 0.78125, 0.5    , 0.375  , 0.9375 , 0.625  , 0.875  ,
        0.90625, 0.40625, 0.65625],
       [0.84375, 0.1875 , 0.4375 , 0.5625 , 0.34375, 0.03125, 0.6875 ,
        0.15625, 0.0625 , 0.75   ]])

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

In [15]:
halton = dr.getHaltonDraws(sampleSize=1,numberOfDraws=10,base=3)
halton
Out[15]:
array([[0.33333333, 0.66666667, 0.11111111, 0.44444444, 0.77777778,
        0.22222222, 0.55555556, 0.88888889, 0.03703704, 0.37037037]])

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

In [16]:
halton = dr.getHaltonDraws(sampleSize=1,numberOfDraws=10,base=3,skip=10)
halton
Out[16]:
array([[0.33333333, 0.66666667, 0.11111111, 0.44444444, 0.77777778,
        0.22222222, 0.55555556, 0.88888889, 0.03703704, 0.37037037]])

Antithetic

Antithetic draws can be generated from any function generating uniform draws

In [17]:
draws = dr.getAntithetic(dr.getUniform,sampleSize=3,numberOfDraws=10)
draws
Out[17]:
array([[0.42212598, 0.47179374, 0.148715  , 0.73860092, 0.04105165,
        0.57787402, 0.52820626, 0.851285  , 0.26139908, 0.95894835],
       [0.51232324, 0.76571874, 0.50131744, 0.1464681 , 0.55853883,
        0.48767676, 0.23428126, 0.49868256, 0.8535319 , 0.44146117],
       [0.41770545, 0.77892061, 0.75202037, 0.04320253, 0.25552776,
        0.58229455, 0.22107939, 0.24797963, 0.95679747, 0.74447224]])

Antithetic MLHS

In [18]:
draws = dr.getAntithetic(dr.getLatinHypercubeDraws,sampleSize=3,numberOfDraws=10)
draws
Out[18]:
array([[0.0780733 , 0.31379531, 0.15866299, 0.0151855 , 0.42181169,
        0.9219267 , 0.68620469, 0.84133701, 0.9848145 , 0.57818831],
       [0.47828697, 0.76997285, 0.24430161, 0.61220294, 0.90674633,
        0.52171303, 0.23002715, 0.75569839, 0.38779706, 0.09325367],
       [0.94935667, 0.56168759, 0.72254999, 0.35044001, 0.80752534,
        0.05064333, 0.43831241, 0.27745001, 0.64955999, 0.19247466]])

Antithetic Halton.

In [19]:
draws = dr.getAntithetic(dr.getHaltonDraws,sampleSize=1,numberOfDraws=10)
draws
Out[19]:
array([[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

In [20]:
def halton(sampleSize,numberOfDraws):
    return dr.getHaltonDraws(numberOfDraws,sampleSize,skip=100)
draws = dr.getAntithetic(halton,sampleSize=3,numberOfDraws=10)
draws
Out[20]:
array([[0.5   , 0.25  , 0.75  , 0.5   , 0.75  , 0.25  ],
       [0.125 , 0.625 , 0.375 , 0.875 , 0.375 , 0.625 ],
       [0.875 , 0.0625, 0.5625, 0.125 , 0.9375, 0.4375],
       [0.3125, 0.8125, 0.1875, 0.6875, 0.1875, 0.8125],
       [0.6875, 0.4375, 0.9375, 0.3125, 0.5625, 0.0625]])

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

In [21]:
draws = dr.getNormalWichuraDraws(sampleSize=3,numberOfDraws=10)
draws
Out[21]:
array([[ 0.46838353,  0.30944461, -0.72705093,  1.47731844,  2.32498619,
        -0.38523969,  2.41753668, -0.89851703,  1.04351099, -0.52799835],
       [ 1.92824184,  0.39265187, -1.69805707,  0.72243106,  0.23153714,
         1.29178853, -0.39065241,  0.75741968, -0.1117659 , -0.25373762],
       [-1.62848672,  0.80968139, -0.47731943, -0.6055917 , -1.04245527,
         2.17223027, -1.20160168,  1.83144122, -1.49289716,  0.7304994 ]])

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

In [22]:
draws = dr.getNormalWichuraDraws(sampleSize=3,numberOfDraws=10,antithetic=True)
draws
Out[22]:
array([[-0.39745017, -0.39736937,  1.09672047, -0.51651723,  0.25290974,
         0.39745017,  0.39736937, -1.09672047,  0.51651723, -0.25290974],
       [-0.3183153 ,  0.48156301,  0.99268525,  0.36495778,  0.60162053,
         0.3183153 , -0.48156301, -0.99268525, -0.36495778, -0.60162053],
       [-0.08619037,  0.78059876, -0.17844824,  0.16459028, -0.03529156,
         0.08619037, -0.78059876,  0.17844824, -0.16459028,  0.03529156]])

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.

In [23]:
myUnif = dr.getLatinHypercubeDraws(sampleSize=3,numberOfDraws=5)
myUnif
Out[23]:
array([[0.73078091, 0.97414181, 0.54479742, 0.0906173 , 0.65406372],
       [0.52895577, 0.22199862, 0.00909876, 0.42051145, 0.87765406],
       [0.15495423, 0.84811817, 0.76065654, 0.33554539, 0.32617466]])
In [24]:
draws = dr.getNormalWichuraDraws(sampleSize=3,numberOfDraws=10,uniformNumbers=myUnif,antithetic=True)
draws
Out[24]:
array([[ 0.61517648,  1.94548717,  0.11252751, -1.33696337,  0.39631515,
        -0.61517648, -1.94548717, -0.11252751,  1.33696337, -0.39631515],
       [ 0.0726452 , -0.76546072, -2.36112653, -0.20058523,  1.16333929,
        -0.0726452 ,  0.76546072,  2.36112653,  0.20058523, -1.16333929],
       [-1.01541415,  1.02839586,  0.70841608, -0.42465146, -0.45050087,
         1.01541415, -1.02839586, -0.70841608,  0.42465146,  0.45050087]])

The same with Halton draws

In [25]:
myUnif = dr.getHaltonDraws(sampleSize=2,numberOfDraws=5,base=3,skip=10)
myUnif
Out[25]:
array([[0.33333333, 0.66666667, 0.11111111, 0.44444444, 0.77777778],
       [0.22222222, 0.55555556, 0.88888889, 0.03703704, 0.37037037]])
In [26]:
draws = dr.getNormalWichuraDraws(numberOfDraws=10,sampleSize=2,uniformNumbers=myUnif,antithetic=True)
draws
Out[26]:
array([[-0.4307273 ,  0.4307273 , -1.22064035, -0.1397103 ,  0.76470967,
         0.4307273 , -0.4307273 ,  1.22064035,  0.1397103 , -0.76470967],
       [-0.76470967,  0.13971031,  1.22064035, -1.78615554, -0.33087257,
         0.76470967, -0.13971031, -1.22064035,  1.78615554,  0.33087257]])
In [ ]: