# Chapter 4: Bayesian Statistics (Working in progress proof of concept)

``````

In [1]:

# These two lines are necessary only if gempy is not installed
import sys, os
sys.path.append("../")

# Importing gempy
import gempy as gp

# Embedding matplotlib figures into the notebooks
%matplotlib inline

# Aux imports
import numpy as np

``````
``````

In [2]:

``````
``````

In [3]:

# Assigning series to formations as well as their order (timewise)
gp.set_data_series(geo_data, {"EarlyGranite_Series": 'EarlyGranite',
"BIF_Series":('SimpleMafic2', 'SimpleBIF'),
"SimpleMafic_Series":'SimpleMafic1'},
order_series = ["EarlyGranite_Series",
"BIF_Series",
"SimpleMafic_Series"], verbose=1)

``````
``````

Out[3]:

text-align: right;
}

text-align: left;
}

.dataframe tbody tr th {
vertical-align: top;
}

EarlyGranite_Series
BIF_Series
SimpleMafic_Series

0
EarlyGranite
SimpleMafic2
SimpleMafic1

1
EarlyGranite
SimpleBIF
SimpleMafic1

``````

Setting uncertainties adding the values to the Dataframe.

``````

In [4]:

geo_data.interfaces['X_std'] = None
geo_data.interfaces['Y_std'] = 0
geo_data.interfaces['Z_std'] = 100

geo_data.foliations['X_std'] = None
geo_data.foliations['Y_std'] = 0
geo_data.foliations['Z_std'] = 0

geo_data.foliations['dip_std'] = 10
geo_data.foliations['azimuth_std'] = 10

``````
``````

Out[4]:

text-align: right;
}

text-align: left;
}

.dataframe tbody tr th {
vertical-align: top;
}

X
Y
Z
azimuth
dip
polarity
formation
series
order_series
G_x
G_y
G_z
X_std
Y_std
Z_std
dip_std
azimuth_std

0
735082.0630
6879102.25
480.551436
276.153239
80.0
1
EarlyGranite
EarlyGranite_Series
1
-0.979134
0.105560
0.173648
None
0
0
10
10

1
715991.2815
6882773.25
505.165864
152.654159
80.0
1
EarlyGranite
EarlyGranite_Series
1
0.452382
-0.874755
0.173648
None
0
0
10
10

2
728767.4065
6878759.25
470.031623
165.980598
80.0
1
EarlyGranite
EarlyGranite_Series
1
0.238570
-0.955474
0.173648
None
0
0
10
10

3
730627.5315
6880472.50
477.402658
120.986348
80.0
1
EarlyGranite
EarlyGranite_Series
1
0.844266
-0.507012
0.173648
None
0
0
10
10

4
732683.4690
6882332.75
481.711952
161.600709
80.0
1
EarlyGranite
EarlyGranite_Series
1
0.310842
-0.934464
0.173648
None
0
0
10
10

``````
``````

In [5]:

# input_data_T = interp_data.interpolator.tg.input_parameters_list()
# select = interp_data.interpolator.pandas_rest_layer_points['formation'] == 'Reservoir'

``````
``````

In [6]:

interp_data = gp.InterpolatorInput(geo_data, compile_theano=False,

``````
``````

I am in the setting
float32
I am here
[2, 2]

``````

Now the generation of the geomodel will be an operation embedded in a larger tree.

``````

In [7]:

import theano
import theano.tensor as T
geomodel = theano.OpFromGraph(interp_data.interpolator.tg.input_parameters_list(),
[interp_data.interpolator.tg.whole_block_model(0)],
on_unused_input='ignore',
)

``````

Because now the GeMpy model is a theano operation and not a theano function, to call it we need to use theano variables (with theano functions we call them with python variables). This is very easy to modify, we just need to use theano shared to convert our python input data into theano variables.

The pymc3 objects are already theano variables (pm.Normal and so on). Now the trick is that using the theano function T.set_subtensor, we can change one deterministic value of the input arrays(the ones printed in the cell above) by a stochastic pymc3 object. Then with the new arrays we just have to call the theano operation and pymc will do the rest

``````

In [8]:

# This is the creation of the model
import pymc3 as pm
theano.config.compute_test_value = 'off'
#theano.config.warn_float64 = 'warn'
model = pm.Model()
with model:
# We create the Stochastic parameters. In this case only the Z position
# of the interfaces
Z_rest = pm.Normal('Z_unc_rest',
interp_data.interpolator.pandas_rest_layer_points['Z'].as_matrix().astype('float32'),
interp_data.interpolator.pandas_rest_layer_points['Z_std'].as_matrix().astype('float32'),
dtype='float32', shape = (66))

Z_ref = pm.Normal('Z_unc_ref', interp_data.interpolator.ref_layer_points[:, 2].astype('float32'),
interp_data.interpolator.ref_layer_points[:, 2].astype('float32'),
dtype='float32', shape = (66))

# We convert a python variable to theano.shared
input_sh = []
for i in interp_data.get_input_data():
input_sh.append(theano.shared(i))

# We add the stochastic value to the correspondant array. rest array is
# a n_points*3 (XYZ) array. We only want to change Z in this case.
input_sh[4] = T.set_subtensor(
input_sh[4][:, 2], Z_ref)

input_sh[5] = T.set_subtensor(
input_sh[5][:, 2], Z_rest)

# With the stochastic parameters we create the geomodel result:
geo_model = pm.Deterministic('GemPy', geomodel(input_sh[0], input_sh[1], input_sh[2],
input_sh[3], input_sh[4], input_sh[5]))

``````
``````

[9, 9, 9]

/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number: 2.4215285776563178e-08
' condition number: {}'.format(rcond), RuntimeWarning)

``````
``````

In [9]:

theano.config.compute_test_value = 'ignore'
# This is the sampling
# BEFORE RUN THIS FOR LONG CHECK IN THE MODULE THEANOGRAF THAT THE FLAG
# THEANO OPTIMIZER IS IN 'fast_run'!!
with model:
# backend = pm.backends.ndarray.NDArray('geomodels')
step = pm.NUTS()
trace = pm.sample(30, tune=10, init=None, step=step, )

``````
``````

0%|          | 0/40 [00:00<?, ?it/s]/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number: 2.4215285776563178e-08
' condition number: {}'.format(rcond), RuntimeWarning)
8%|▊         | 3/40 [00:04<01:18,  2.13s/it]/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number: 3.6804916447863434e-08
' condition number: {}'.format(rcond), RuntimeWarning)
10%|█         | 4/40 [00:05<01:00,  1.69s/it]/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number: 4.177091028623181e-08
' condition number: {}'.format(rcond), RuntimeWarning)
12%|█▎        | 5/40 [00:06<00:48,  1.38s/it]/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number: 2.817358790707658e-08
' condition number: {}'.format(rcond), RuntimeWarning)
15%|█▌        | 6/40 [00:06<00:39,  1.16s/it]/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number: 5.478779740997197e-08
' condition number: {}'.format(rcond), RuntimeWarning)
18%|█▊        | 7/40 [00:07<00:33,  1.01s/it]/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number: 4.4421348377454706e-08
' condition number: {}'.format(rcond), RuntimeWarning)
20%|██        | 8/40 [00:08<00:28,  1.11it/s]/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number: 5.4531707149862996e-08
' condition number: {}'.format(rcond), RuntimeWarning)
28%|██▊       | 11/40 [00:10<00:21,  1.36it/s]/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number: 4.235039341438096e-08
' condition number: {}'.format(rcond), RuntimeWarning)
40%|████      | 16/40 [00:13<00:15,  1.51it/s]/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number: 5.433537708654512e-08
' condition number: {}'.format(rcond), RuntimeWarning)
42%|████▎     | 17/40 [00:14<00:15,  1.52it/s]/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number: 4.1121396066046145e-08
' condition number: {}'.format(rcond), RuntimeWarning)
45%|████▌     | 18/40 [00:14<00:14,  1.52it/s]/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number: 4.794599561819268e-08
' condition number: {}'.format(rcond), RuntimeWarning)
85%|████████▌ | 34/40 [00:25<00:03,  1.53it/s]/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number: 3.8119761569532784e-08
' condition number: {}'.format(rcond), RuntimeWarning)
92%|█████████▎| 37/40 [00:27<00:01,  1.51it/s]/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number: 4.07877926988931e-08
' condition number: {}'.format(rcond), RuntimeWarning)
98%|█████████▊| 39/40 [00:28<00:00,  1.47it/s]/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number: 3.935079107009187e-08
' condition number: {}'.format(rcond), RuntimeWarning)
100%|██████████| 40/40 [00:29<00:00,  1.47it/s]/home/miguel/anaconda3/lib/python3.6/site-packages/pymc3/step_methods/hmc/nuts.py:418: UserWarning: Chain 0 contains only 30 samples.
% (self._chain_id, n))
/home/miguel/anaconda3/lib/python3.6/site-packages/pymc3/step_methods/hmc/nuts.py:448: UserWarning: Chain 0 reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.
'reparameterize.' % self._chain_id)

``````
``````

In [10]:

import matplotlib.pyplot as plt
for i in range(100):
gp.plot_section(geo_data, trace.get_values('GemPy')[i][-1, 0, :], 13,
direction='y', plot_data=False)
plt.show()

``````
``````

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-10-bb562bc0d205> in <module>()
1 import matplotlib.pyplot as plt
2 for i in range(100):
----> 3     gp.plot_section(geo_data, trace.get_values('GemPy')[i][-1, 0, :], 13,
4                        direction='y', plot_data=False)
5     plt.show()

IndexError: index 30 is out of bounds for axis 0 with size 30

``````
``````

In [ ]:

``````