Analysis of uncertainty propagation with Monte Carlo sampling



In [1]:

    
from IPython.core.display import HTML
css_file = 'pynoddy.css'
HTML(open(css_file, "r").read())









    Out[1]:



In [3]:

    
%matplotlib inline



In [36]:

    
# here the usual imports. If any of the imports fails, make sure that pynoddy is installed
# properly, ideally with 'python setup.py develop' or 'python setup.py install'
import sys, os
import matplotlib.pyplot as plt
import numpy as np
# adjust some settings for matplotlib
from matplotlib import rcParams
# print rcParams
rcParams['font.size'] = 15
# determine path of repository to set paths corretly below
repo_path = os.path.realpath('../..')
import pynoddy.history
import pynoddy.experiment
rcParams.update({'font.size': 20})

Using the Monte-Carlo experiment class

The methods for the analysis of uncertainty propagation with Monte-Carlo sampling are integrated in a Python experiment class. As before, we can therefore load the model and a parameter file which defines the parameter statistics for the sampling step:



In [162]:

    
reload(pynoddy.history)
reload(pynoddy.output)
reload(pynoddy.experiment)
reload(pynoddy)
# the model itself is now part of the repository, in the examples directory:
history_file = os.path.join(repo_path, "examples/GBasin_Ve1_V4_b.his")
gipps_mc = pynoddy.experiment.MonteCarlo(history = history_file)









    



 STRATIGRAPHY
 FOLD
 UNCONFORMITY
 FAULT
 FAULT
 UNCONFORMITY
 FAULT
 FAULT
 UNCONFORMITY
 FAULT
 FAULT
 UNCONFORMITY
 TILT
 FOLD



In [163]:

    
gipps_mc.load_parameter_file(os.path.join(repo_path, "examples/gipps_params_2.csv"))
gipps_mc.freeze()

We now define the number of samples to be generated:



In [164]:

    
gipps_mc.set_n_draws(100)
probs = gipps_mc.perform_sampling('', output_type = 'section', resolution = 20)



In [165]:

    
probs.keys()









    Out[165]:





[1.0,
 2.0,
 3.0,
 4.0,
 5.0,
 7.0,
 10.0,
 11.0,
 12.0,
 13.0,
 14.0,
 15.0,
 16.0,
 17.0,
 18.0,
 19.0,
 20.0]



In [166]:

    
gipps_mc.plot_section(direction = 'x', ve = 5.0)



In [167]:

    
gipps_mc.nx = np.shape(probs[7.0])[0]
gipps_mc.nz = np.shape(probs[7.0])[2]
gipps_mc.plot_section(direction = 'x', data = probs[7.0], ve = 5., cmap = 'RdBu_r')

Furthermore, we set the random seed to enable a completely reproducible experiment:



In [52]:

    
gipps_mc.set_random_seed(180415)



In [49]:

    
fig = plt.figure(figsize = (12,6))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)
gipps_mc.plot_section(ax = ax1, direction = 'x', model_type = "base", 
                           colorbar = False, title = "", ve = 5.)
gipps_mc.plot_section(ax = ax2, direction = 'x', colorbar = False, 
                           title = "", ve = 5.)