# -*- coding: utf-8 -*-
"""
Created on Fri Apr 10 12:02:35 2020

@author: neil
"""

# 1st in a series, First we discuss random numbers, and PDFs in numpy
# in preparation to input PDFs to a Monte Carlo approach.

# This program illustrates sampling from a 'random' distribution to produce
# samples of a random variable.  If you take a lot of samples, and plot the
# histogram of the result, you can sketch the underlying distribution.

# We can view this as finding values of rock conductivity, if we only know the PDFs
# of the conductivities, and want to get representative samples.

import matplotlib.pyplot as plt
import numpy as np

n = 50000                       # number of random samples to take from distribution

fig = plt.figure()              # set up a plot window
fig.canvas.set_window_title ("PDFs of the rock properties")
ax1 = fig.add_subplot(2,1,1)
ax1.grid(True)

# for our current problem we can consider a to represent clays, and b to represent sands
# However, in general the choice of 'distributions' is a deeply difficult part of a problem.
# most people just assume a 'Normal' distribution, but this is often incorrect
# Many sampling problems in the GeoSciences actually have a LogNormal distribution.
a=np.random.lognormal(0.3,.5,size=n) # sample the lognormal distribution
# There are many other possible distributions, but if we don't know much about
# the actual distribution, we can use something simple such as a triangular distribution
b=np.random.triangular(1,4,5,size=n) # sample a triangular distribution (left middle and right points of triangle)

ax1.hist(a,bins=200)        # this calls np.histogram to calculate the bins in a bar plot
ax1.set_xlim([0,6])         # the lognormal has a LONG positive tail


ax1.hist(b,bins=100,histtype='step')  # plot the histogram as an 'transparent' bar plot
ax1.set_ylabel("number in each conductivity bin")
#ax1.set_xlabel("Values of Conductivity")
ax1.set_title("Raw conductivity, random values of 2 rock types")

# to turn the histogram of the distributions into PDFs we make the area under the curves equal 1
# this allows us to interpret the bin y values as probabilities of the conductivity being within the bin.
# 'PDF plots have as the vertical axis not probability, but probability per unit on the x axis.  
# It only makes sense to ask for the probability of a specified range ot the x axis events'
ax2 = fig.add_subplot(2,1,2) # plot the resulting PDFs
ax2.grid(True)
ax2.hist(a,bins=200,density='True')     # this calls np.histogram to calculate the bins, and 'Normalize'
ax2.set_xlim([0,6])                     # the lognormal has a LONG positive tail
ax2.hist(b,bins=100,density='True',histtype='step')
# the triangular distribution is plotted as only a line so we can see the lognormal
ax2.set_ylabel("Probability (density)")
ax2.set_xlabel("Values of Conductivity")
ax2.set_title("PDF, or Normalized distribution, of values of 2 rock types")


plt.show()