This package provides a library to perform basic monte carlo estimation of property distributions in native python. The package was originaly intended to estimate oil and gas reserves in multi-regional domains from petrophysical datasets, but can be applied to a variety of problems such as the example problems shown here.
%matplotlib inline
from montepetro.models import Model
from montepetro.regions import Region
from montepetro.properties import RandomProperty
from montepetro.seed_generators import SeedGenerator
import numpy
MontePetro allows to create single and multi-region monte carlo estimates of parameter distributions. The basic container for all regions is provided by the model class.
seed = 300
model = Model("A simple Model", seed)
The model contains regions as well as properties. It also handles the generation of seed values for all the regional property distributions. See the Model class documentation for more information.
Regions in MontePetro contain all of the properties of the model and have their own distributional properties e.g. two different types of bone structure in a larger bone could be modeled as two regions in a model.
We can perform calculations on the regions themselves or we can pass them to a model structure and let the model do all the hard work for us!
region_a = Region(name="Inner Bone")
region_b = Region(name="Outer Bone")
#For simplicity let's add these to our model
model.add_region(region_a)
model.add_region(region_b)
MontePetro already comes with a number of properties types that we can define for our model. RandomProperties are computed contain sets of values that we can sample from probability distributions. RegionalProperties are properties that we want to calculate based on the defined properties of our model.
For our bone example me way want to estimate the ensemble distribution of the density of the bone. We define a simple model for the ensemble density of the bone to be:
$$\rho_{bone} = \rho_{bone}(1-\phi)+\rho_{fluid}\phi$$The following code goes into detail on how we create the probability distributions and add them to our model.
#Define number of samples we want to take
n = 100000
#We pass the RandomProperty our functions from which we sample the probability distributions
#Here we make use of the numpy libraries built in random number generators
porosity = RandomProperty(name="Porosity", n=n, \
random_number_function=numpy.random.uniform)
density_calcite = RandomProperty(name="Rho_Calc", n=n, \
random_number_function=numpy.random.triangular)
density_bone_fluid = RandomProperty(name="Rho_Fluid", n=n, \
random_number_function=numpy.random.triangular)
config = {"Inner Bone": {"Porosity": {"low": 0.0, "high": 0.2},
"Rho_Calc": {"left": 2.54, "right": 2.7, "mode": 2.58},
"Rho_Fluid": {"left": 1.01, "right": 1.5, "mode": 1.2}},
"Outer Bone": {"Porosity": {"low": 0.1, "high": 0.2},
"Rho_Calc": {"left": 2.50, "right": 3.0, "mode": 2.8},
"Rho_Fluid": {"left": 1.1, "right": 1.3, "mode": 1.2}}}
#Let's add these to the model.
model.add_property(porosity)
model.add_property(density_calcite)
model.add_property(density_bone_fluid)
#Some Model container magic! We add all these properties to the regions.
model.add_defined_properties_to_regions()
#We pass the model our configuration and run the model
#This will generate all the sampled distributions for each region
model.run(config)
We will now plot the ensemble density $\rho_{ensemble}$ as a histogram and evaluate the results.
import matplotlib.pyplot as plot
densities = []
for region_name, region in model.regions.iteritems():
porosity = region.properties["Porosity"].values
rho_calc = region.properties["Rho_Calc"].values
rho_bone_fluid = region.properties["Rho_Fluid"].values
ensemble_density = rho_calc*(1-porosity)+rho_bone_fluid*porosity
densities.append(ensemble_density)
#And we can plot these values of course for our example.
#plot.hist(ensemble_density, bins=50)
total_density = numpy.add(densities[0], densities[1])
plot.hist(total_density, bins=500)
plot.xlabel(u"$Bone \ Density$", fontsize=14)
plot.ylabel(u"$Frequency$",fontsize=14)
plot.show()
This concludes this example. If you find MontePetro useful in your work or research please consider citing it in your references.
To cite MontePetro in publications use:
"MontePetro Development Team (2015). MontePetro: Python library for probabilistic reserve estimates URL http://www.github.com/lukasmosser/MontePetro"