# -*- coding: utf-8 -*-
"""
Created 2022
@author: neil
"""

# 0th in a Monte Carlo series, First we discuss random numbers, and PDF, CDFs, and PMFs in numpy
# in preparation to input PDFs or PMFs 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 find a sketch of the underlying distribution.
# In addition, you can convert a histogram of the distribution of the discrete random
# numbers to make a 'probability mass function' (PMF) by normalizing.
# Normalizing means making the sum of the distribution equal to one.

# Although CDFs and PDFs are usually more associated with 'continuous' random variables
# we can also illustrates making a cummulative density function (CDF), which is the
# integral of the PMF (or PDF), and finally a probability density function (PDF)
# which is the derivative of the (smoothed) CDF.

import matplotlib.pyplot as plt
import numpy as np

n   = 100000                       # number of random samples to take from distribution
bins= 200                          # number of bins for histograms for plotting

# For our current problem we consider the distribution of thermal conductivities of clays.
# To illustrate we consider the conductivities to have a 'lognormal distribution' about a center.
# 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, especially in
# geoscience problems. We use numpy.random to get a 'sample' from the lognormal distribution.

fig = plt.figure()                    # set up a plot window, and plot the lognormal sample
ax1 = fig.add_subplot(2,2,1)
ax1.grid(True)
A = np.random.lognormal(0.3,.5,size=n)# sample the lognormal distribution, A is a vector length 'n'
ax1.hist(A,bins=bins)                  # this calls np.histogram to count the values into the 200 bins
ax1.set_xlim([0,6])                   # the lognormal has a LONG positive tail, so truncate at '6'
ax1.set_ylabel("number in each conductivity bin")
ax1.set_title("Raw bin count for lognormal conductivity")

# to turn the histogram of the distributions into PMFs we 'Normalize' or divide by numbers in bins
# PMF 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 or spread of the x axis events
ax2 = fig.add_subplot(2,2,3)
ax2.grid(True)
ax2.hist(A,bins=bins,density='True')     # this calls np.histogram to calculate the bins
ax2.set_xlim([0,6])                     # the lognormal has a LONG positive tail
ax2.set_ylabel("Probability (mass density)")
ax2.set_xlabel("Values of Conductivity")
ax2.set_title("PMF, or Normalized distribution")

# A CDF (cummulative Density Function) describes the cummulative probability that the
# random number will be found to be less than a value on the x axis.
ax3 = fig.add_subplot(2,2,2)
ax3.grid(True)
B = np.histogram(A,bins=bins) # uses the same random numbers (A) to make just the bin counts, no plot 
b = np.array(B[0])           # B is 2 lists, one of the numbers, the other the bin locations
xb= np.array(B[1])           # so b is the bin counts, and xb are their locations
# actually the middle of each bin 'i' is at (xb[i]+xb[i+1])/2
S = np.zeros(len(b))
S[0] = b[0]
xm = np.zeros(len(b))
xm[0] = (xb[0]+xb[1])/2.
for i in range(1,len(b)):    # sum the number of samples that are less_equal than each bin
    S[i] = S[i-1]+b[i]       # S contains the number of samples less_equal to each bin 'i'
    xm[i]= (xb[i]+xb[i+1])/2.# also use this loop to make center locations for each bin
C = S/n                      # by dividing by the total number of samples CDF[-1] will be 1
ax3.plot(xm,C)               # plot the NORMALIZED sum 'C', which is the CDF, going from 0 to 1
ax3.set_xlim([0,6])          # the lognormal has a LONG positive tail, so truncate in plot
ax3.set_ylabel("Probability of value being less than:")
ax3.set_title("CDF, or cummulative probability")

# if we take the derivative of the CDF we get the PDF, although in this case it is only an
# estimate since we are deriving it from discrete random numbers, not from the underlying distribution
ax4 = fig.add_subplot(2,2,4)
ax4.grid(True)

# smooth C before taking derivative (useful idea, differentiating amplifies noise!)
# A simple 'box car' smoothing
CDF = np.zeros_like(C)
CDF[0] = C[0]
CDF[-1]= C[-1]
for i in range(1,len(C)-1):
    CDF[i] = (C[i-1] + C[i] + C[i+1])/3
# the CDF has now been somewhat smoothed, and we can differentiate, using np.diff

PDF = np.diff(CDF)                      # np.diff returns 1 less value than is sent to it
PDF = PDF/(xb[1]-xb[0])                 # to make the difference array a derivative, divide by dx
ax4.plot(xb[1:-1],PDF,'r')              # minor point, the PDF locations are actually at xb, not xm
ax4.plot([xb[0],xb[1]],[0,PDF[0]],'r')  # another minor detail, since PDF doesn't include the zero
ax4.set_xlim([0,6])                     # the lognormal has a LONG positive tail
ax4.set_ylabel("Probability density")
ax4.set_title("PDF, derived from sampling")
ax4.set_xlabel("Values of Conductivity")
# as a final check on all our calculations, the integral of the PDF should = 1
print("Integral of PDF (should be 1) ",np.sum(PDF)*(xb[1]-xb[0]))

plt.show()