<p align="center">
    

</p>

## Data Analytics 

### Parametric Distributions  in Python 


#### Michael Pyrcz, Associate Professor, University of Texas at Austin 

##### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)


### Data Analytics: Parametric Distributions

Here's a demonstration of making and general use of parametric distributions in Python. This demonstration is part of the resources that I include for my courses in Spatial / Subsurface Data Analytics at the Cockrell School of Engineering at the University of Texas at Austin.  

#### Parametric Distributions

We will cover the following distributions:

* Uniform
* Triangular
* Gaussian
* Log Normal

We will demonstrate:

* distribution parameters 
* forward and inverse operators
* summary statistics

I have a lecture on these parametric distributions available on [YouTube](https://www.youtube.com/watch?v=U7fGsqCLPHU&t=1687s).   

#### Getting Started

Here's the steps to get setup in Python with the GeostatsPy package:

1. Install Anaconda 3 on your machine (https://www.anaconda.com/download/). 
2. From Anaconda Navigator (within Anaconda3 group), go to the environment tab, click on base (root) green arrow and open a terminal. 
3. In the terminal type: pip install geostatspy. 
4. Open Jupyter and in the top block get started by copy and pasting the code block below from this Jupyter Notebook to start using the geostatspy functionality. 

You will need to copy the data file to your working directory.  They are available here:

* Tabular data - unconv_MV_v4.csv at https://git.io/fhHLT.

#### Importing Packages

We will need some standard packages. These should have been installed with Anaconda 3.


```python
import numpy as np                        # ndarrys for gridded data
import pandas as pd                       # DataFrames for tabular data
import os                                 # set working directory, run executables
import matplotlib.pyplot as plt           # for plotting
from scipy import stats                   # summary statistics
import math                               # trigonometry etc.
import scipy.signal as signal             # kernel for moving window calculation
import random                             # for randon numbers
import seaborn as sns                     # for matrix scatter plots
from scipy import linalg                  # for linear regression
from sklearn import preprocessing
```

#### Set the Working Directory

I always like to do this so I don't lose files and to simplify subsequent read and writes (avoid including the full address each time). 


```python
os.chdir("c:/PGE383")                     # set the working directory
```

### Uniform Distribution

Let's start with the most simple distribution.

* by default a random number is uniform distributed

* this ensures that enough random samples (Monte Carlo simulations) will reproduce the distribution

\begin{equation}
x_{\alpha}^{s} = F^{-1}_x(p_{\alpha}), \quad X^{s} \sim F_X
\end{equation}

#### Random Samples

Let's demonstrate the use of the command:

```python
uniform.rvs(size=n, loc = low, scale = interval, random_state = seed)
```

Where:

* size is the number of samples

* loc is the minimum value

* scale is the range, maximum value minus the minimum value

* random_state is the random number seed

We will observe the convergence of the samples to a uniform distribution as the number of samples becomes large.

We will make a compact set of code by looping over all the cases of number of samples

* we store the number of samples cases in the list called ns

* we store the samples as a list of lists, called X_uniform



```python
from scipy.stats import uniform

low = 0.05; interval = 0.20; ns = [1e1,1e2,1e3,1e4,1e5,1e6]; X_uniform = []

index = 0
for n in ns:
    X_uniform.append(uniform.rvs(size=int(ns[index]), loc = low, scale = interval).tolist())
    plt.subplot(2,3,index+1)
    GSLIB.hist_st(X_uniform[index],loc,loc+interval,log=False,cumul = False,bins=20,weights = None,xlabel='Values',title='Distribution, N = ' + str(int(ns[index])))
    index = index + 1
    
plt.subplots_adjust(left=0.0, bottom=0.0, right=2.3, top=1.6, wspace=0.2, hspace=0.3)
```

We can observe that by drawing more Monte Carlo simulations, we more closely approximate the original uniform parametric distribution.

#### Forward Distribution

Let's demonstrate the forward operator. We can take any value and calculate the associated:

* density (probability density function)
* cumulative probability

The transform for the probability density function is:

\begin{equation}
p = f_x(x)
\end{equation}

where $f_x$ is the PDF and $p$ is the density for value, $x$.

and for the cumulative distribution function is:

\begin{equation}
P = F_x(x)
\end{equation}

where $F_x$ is the CDF and $P$ is the cumulative probability for value, $x$.


```python
x_values = np.linspace(0.0,0.3,100)
p_values = uniform.pdf(x_values, loc = low, scale = interval)
P_values = uniform.cdf(x_values, loc = low, scale = interval)

plt.subplot(1,2,1)
plt.plot(x_values, p_values,'r-', lw=5, alpha=0.3, label='uniform PDF'); plt.title('Uniform PDF'); plt.xlabel('Values'); plt.ylabel('Density')

plt.subplot(1,2,2)
plt.plot(x_values, P_values,'r-', lw=5, alpha=0.3, label='uniform CDF'); plt.title('Uniform CDF'); plt.xlabel('Values'); plt.ylabel('Cumulative Probability')

plt.subplots_adjust(left=0.0, bottom=0.0, right=1.8, top=0.8, wspace=0.2, hspace=0.3)
```

#### Inverse Distribution

Let's know demonstrate the reverse operator for the uniform distribution:

\begin{equation}
X = F^{-1}_X(P)
\end{equation}


```python
p_values = np.linspace(0.01,0.99,100)
x_values = uniform.ppf(p_values, loc = low, scale = interval)
plt.plot(x_values, p_values,'r-', lw=5, alpha=0.3, label='uniform pdf')
```

#### Summary Statistics

We also have a couple of convience member functions to return the statistics from the parametric distribution:

* mean
* median
* mode
* variance
* standard deviation

Let's demonstrate a few of these methods.

```python
uniform.stats(loc = low, scale = interval, moments = 'mvsk')
```

returns a tuple with the mean, variance, skew and kurtosis (centered 1st, 2nd, 3rd and 4th moments)


```python
print('Stats: mean, variance, skew and kurtosis = ' + str(uniform.stats(loc = low, scale = interval, moments = 'mvsk')))
```

    Stats: mean, variance, skew and kurtosis = (array(0.15), array(0.00333333), array(0.), array(-1.2))


We can confirm this by calculating the centered variance (regular variance) with this member function:

```python
uniform.var(loc = low, scale = interval)
```


```python
print('The variance is ' + str(round(uniform.var(loc = low, scale = interval),4)) + '.')
```

    The variance is 0.0033.


We can also directly calculate the:

* standard deviation - std
* mean - mean
* median - median

We can also calculate order of a non-centered moment. The moment method allows us to calculate an non-centered moment of any order. Try this out.


```python
m_order = 4
print('The ' + str(m_order) + 'th order non-centered moment is ' + str(uniform.moment(n = m_order, loc = low, scale = interval)))
```

    The 4th order non-centered moment is 0.0009762500000000001


#### Symmetric Interval

We can also get the symmetric interval (e.g. prediction or confidence intervals) for any alpha level. 

* Note the program mislabels the value as alpha, it is actually the significance level (1 - alpha)


```python
level = 0.95
print('The interval at alpha level ' + str(round(1-level,3)) + ' is ' + str(uniform.interval(alpha = alpha,loc = low,scale = interval)))
```

    The interval at alpha level 0.05 is (0.05500000000000001, 0.245)


#### Triangular Distribution

The great thing about parametric distributions is that the above member functions are the same!

* we can plug and play other parametric distributions and repeat the above.

This time we will make it much more compact!  

* we will import the triangular distribution as my_dist and call the same functions as before
* we need a new parameter, the distribution mode (c parameter)


```python
from scipy.stats import triang as my_dist                                    # import traingular dist as my_dist
dist_type = 'Triangular'                                                     # give the name of the distribution for labels
low = 0.05; mode = 0.20; c = 0.10                                            # given the distribution parameters

x_values = np.linspace(0.0,0.3,100)                                          # get an array of x values
p_values = my_dist.pdf(x_values, loc = low, c = mode, scale = interval)      # calculate density for each x value
P_values = my_dist.cdf(x_values, loc = low, c = mode, scale = interval)      # calculate cumulative probablity for each x value

plt.subplot(1,3,1)                                                           # plot the resulting PDF
plt.plot(x_values, p_values,'r-', lw=5, alpha=0.3); plt.title('Sampling ' + str(dist_type) + ' PDF'); plt.xlabel('Values'); plt.ylabel('Density')

plt.subplot(1,3,2)                                                           # plot the resulting CDF
plt.plot(x_values, P_values,'r-', lw=5, alpha=0.3); plt.title('Sampling ' + str(dist_type) + ' CDF'); plt.xlabel('Values'); plt.ylabel('Cumulative Probability')

p_values = np.linspace(0.00001,0.99999,100)                                 # get an array of p-values
x_values = my_dist.ppf(p_values, loc = low, c = mode, scale = interval)      # apply inverse to get x values from p-values
plt.subplot(1,3,3)
plt.plot(x_values, p_values,'r-', lw=5, alpha=0.3, label='uniform pdf')

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.8, top=0.8, wspace=0.2, hspace=0.3); plt.title('Sampling Inverse ' + str(dist_type) + ' CDF'); plt.xlabel('Values'); plt.ylabel('Cumulative Probability')

print('The mean is ' + str(round(uniform.mean(loc = low, scale = interval),4)) + '.') # calculate stats and symmetric interval
print('The variance is ' + str(round(uniform.var(loc = low, scale = interval),4)) + '.')
print('The interval at alpha level ' + str(round(1-level,3)) + ' is ' + str(uniform.interval(alpha = alpha,loc = low,scale = interval)))
```

#### Gaussian Distribution

Let's now use the Gaussian parametric distribution.

* we will need the parameters mean and the variance

We will apply the forward and reverse operations and calculate the summary statistics.



```python
from scipy.stats import norm as my_dist                                      # import traingular dist as my_dist
dist_type = 'Gaussian'                                                       # give the name of the distribution for labels
mean = 0.15; stdev = 0.05                                                    # given the distribution parameters

x_values = np.linspace(0.0,0.3,100)                                          # get an array of x values
p_values = my_dist.pdf(x_values, loc = mean, scale = stdev)                  # calculate density for each x value
P_values = my_dist.cdf(x_values, loc = mean, scale = stdev)                  # calculate cumulative probablity for each x value

plt.subplot(1,3,1)                                                           # plot the resulting PDF
plt.plot(x_values, p_values,'r-', lw=5, alpha=0.3); plt.title('Sampling ' + str(dist_type) + ' PDF'); plt.xlabel('Values'); plt.ylabel('Density')

plt.subplot(1,3,2)                                                           # plot the resulting CDF
plt.plot(x_values, P_values,'r-', lw=5, alpha=0.3); plt.title('Sampling ' + str(dist_type) + ' CDF'); plt.xlabel('Values'); plt.ylabel('Cumulative Probability')

p_values = np.linspace(0.00001,0.99999,100)                                  # get an array of p-values
x_values = my_dist.ppf(p_values, loc = mean, scale = stdev)                  # apply inverse to get x values from p-values
plt.subplot(1,3,3)
plt.plot(x_values, p_values,'r-', lw=5, alpha=0.3, label='uniform pdf')

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.8, top=0.8, wspace=0.2, hspace=0.3); plt.title('Sampling Inverse ' + str(dist_type) + ' CDF'); plt.xlabel('Values'); plt.ylabel('Cumulative Probability')

print('The mean is ' + str(round(my_dist.mean(loc = mean, scale = stdev),4)) + '.') # calculate stats and symmetric interval
print('The variance is ' + str(round(my_dist.var(loc = mean, scale = stdev),4)) + '.')
print('The interval at alpha level ' + str(round(1-level,3)) + ' is ' + str(my_dist.interval(alpha = alpha,loc = mean,scale = stdev)))
```

#### Log Normal Distribution

Now let's check out the log normal distribution.  

* We need the parameters $\mu$ and $\sigma$


```python
from scipy.stats import lognorm as my_dist                                   # import traingular dist as my_dist
dist_type = 'Log Normal'                                                     # give the name of the distribution for labels
mu = np.log(0.10); sigma = 0.2                                               # given the distribution parameters

x_values = np.linspace(0.0,0.3,100)                                          # get an array of x values
p_values = my_dist.pdf(x_values, s = sigma, scale = np.exp(mu))              # calculate density for each x value
P_values = my_dist.cdf(x_values, s = sigma, scale = np.exp(mu))              # calculate cumulative probablity for each x value

plt.subplot(1,3,1)                                                           # plot the resulting PDF
plt.plot(x_values, p_values,'r-', lw=5, alpha=0.3); plt.title('Sampling ' + str(dist_type) + ' PDF'); plt.xlabel('Values'); plt.ylabel('Density')

plt.subplot(1,3,2)                                                           # plot the resulting CDF
plt.plot(x_values, P_values,'r-', lw=5, alpha=0.3); plt.title('Sampling ' + str(dist_type) + ' CDF'); plt.xlabel('Values'); plt.ylabel('Cumulative Probability')

p_values = np.linspace(0.00001,0.99999,100)                                        # get an array of p-values
x_values = my_dist.ppf(p_values,  s = sigma, scale = np.exp(mu))             # apply inverse to get x values from p-values
plt.subplot(1,3,3)
plt.plot(x_values, p_values,'r-', lw=5, alpha=0.3, label='uniform pdf')

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.8, top=0.8, wspace=0.2, hspace=0.3); plt.title('Sampling Inverse ' + str(dist_type) + ' CDF'); plt.xlabel('Values'); plt.ylabel('Cumulative Probability')

#print('The mean is ' + str(round(my_dist.mean(loc = mean, scale = stdev),4)) + '.') # calculate stats and symmetric interval
#print('The variance is ' + str(round(my_dist.var(loc = mean, scale = stdev),4)) + '.')
#print('The interval at alpha level ' + str(round(1-level,3)) + ' is ' + str(my_dist.interval(alpha = alpha,loc = mean,scale = stdev)))
```

There are many other parametric distributions that we could have included.  Also we could have demonstrated the distribution fitting. 

#### Comments

This was a basic demonstration of working with parametric distributions. 

I have other demonstrations on the basics of working with DataFrames, ndarrays, univariate statistics, plotting data, declustering, data transformations, trend modeling and many other workflows available at [Python Demos](https://github.com/GeostatsGuy/PythonNumericalDemos) and a Python package for data analytics and geostatistics at [GeostatsPy](https://github.com/GeostatsGuy/GeostatsPy). 
  
I hope this was helpful,

*Michael*

#### The Author:

### Michael Pyrcz, Associate Professor, University of Texas at Austin 
*Novel Data Analytics, Geostatistics and Machine Learning Subsurface Solutions*

With over 17 years of experience in subsurface consulting, research and development, Michael has returned to academia driven by his passion for teaching and enthusiasm for enhancing engineers' and geoscientists' impact in subsurface resource development. 

For more about Michael check out these links:

#### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)

#### Want to Work Together?

I hope this content is helpful to those that want to learn more about subsurface modeling, data analytics and machine learning. Students and working professionals are welcome to participate.

* Want to invite me to visit your company for training, mentoring, project review, workflow design and / or consulting? I'd be happy to drop by and work with you! 

* Interested in partnering, supporting my graduate student research or my Subsurface Data Analytics and Machine Learning consortium (co-PIs including Profs. Foster, Torres-Verdin and van Oort)? My research combines data analytics, stochastic modeling and machine learning theory with practice to develop novel methods and workflows to add value. We are solving challenging subsurface problems!

* I can be reached at mpyrcz@austin.utexas.edu.

I'm always happy to discuss,

*Michael*

Michael Pyrcz, Ph.D., P.Eng. Associate Professor The Hildebrand Department of Petroleum and Geosystems Engineering, Bureau of Economic Geology, The Jackson School of Geosciences, The University of Texas at Austin

#### More Resources Available at: [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)



```python

```

1f5e3d1f9e0b6444733670755a29258bec3630fb
ipynb
Jupyter Notebook
PythonDataBasics_ParametricDistributions.ipynb
ushareng/PythonNumericalDemos
a78fba8234760c6449de4a2e5c66348ad9f6e61c
MIT
2020-09-29T19:33:08.000Z
2021-08-08T15:22:21.000Z
PythonDataBasics_ParametricDistributions.ipynb
macca69/PythonNumericalDemos
5c31946198d5735a6262eb63dc75ba11534b49e0
MIT
PythonDataBasics_ParametricDistributions.ipynb
macca69/PythonNumericalDemos
5c31946198d5735a6262eb63dc75ba11534b49e0
MIT
2020-09-26T09:27:16.000Z
2020-11-10T16:40:33.000Z
Qwen/Qwen-72B
1. YES
2. YES