PlotDistributions

PlotDistributions.py

# ------------------------------------------------------------------------
# The following Python code is implemented by Professor Terje Haukaas at
# the University of British Columbia in Vancouver, Canada. It is made
# freely available online at terje.civil.ubc.ca together with notes,
# examples, and additional Python code. Please be cautious when using
# this code; it may contain bugs and comes without warranty of any kind.
# ------------------------------------------------------------------------


# Import useful libraries
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, lognorm, uniform, gamma
from scipy.special import gamma as gammafunction

# Set mean and standard deviation of the random variable
meanX = 3.0
stdvX = 1.0


# Create plot
fig, ax = plt.subplots()


# Lay out the x-axis for plots
x = np.linspace(0.01, 4.0*meanX, 50)


# Normal PDF
normal_pdf_expression = 1.0 / np.sqrt(2 * np.pi  * stdvX**2) * np.exp(-0.5*((x - meanX) / stdvX)**2)
normal_pdf_scipy = norm.pdf(x, meanX, stdvX)
ax.plot(x, normal_pdf_expression, 'ko', label='Normal expression')
ax.plot(x, normal_pdf_scipy, 'k', label='Normal by Scipy')


# Lognormal PDF (Y is the underlying Normal variable)
meanY = np.log(meanX) - 0.5 * np.log(1.0 + (stdvX/meanX)**2)
stdvY = np.sqrt(np.log((stdvX/meanX)**2 + 1.0))
lognormal_pdf_expression = 1.0 / (x * np.sqrt(2 * np.pi  * stdvY**2)) * np.exp(-0.5*((np.log(x) - meanY)/stdvY)**2)
lognormal_pdf_scipy = lognorm.pdf(x, stdvY, 0.0, np.exp(meanY))
ax.plot(x, lognormal_pdf_expression, 'ro', label='Lognormal expression')
ax.plot(x, lognormal_pdf_scipy, 'r', label='Lognormal by Scipy')


# Uniform PDF
a = meanX - np.sqrt(3.0) * stdvX
b = meanX + np.sqrt(3.0) * stdvX
uniform_pdf_expression = np.zeros(len(x))
for i in range(len(x)):
    if x[i] > a and x[i] < b:
        uniform_pdf_expression[i] = 1.0 / (b-a)
uniform_pdf_scipy = uniform.pdf(x, a, b-a)
ax.plot(x, uniform_pdf_expression, 'go', label='Uniform expression')
ax.plot(x, uniform_pdf_scipy, 'g', label='Uniform by Scipy')


# Gamma PDF
a = (meanX / stdvX)**2.0
b = stdvX**2.0 / meanX
k = a
nu = 1/b
gamma_pdf_expression = 1.0 / gammafunction(k) * nu * (nu * x)**(k-1.0) * np.exp(-nu * x)
gamma_pdf_scipy = gamma.pdf(x, a, 0.0, b)
ax.plot(x, gamma_pdf_expression, 'bo', label='Gamma expression')
ax.plot(x, gamma_pdf_scipy, 'b', label='Gamma by Scipy')


# Display the plot
legend = ax.legend(loc='upper right')
plt.show()