There are so many different continuous distributions available in Stan that I often find it hard to see the forest for the trees--- for this reason I've put together a notebook with a cell for every possible continuous distribution available in Stan.
Each function is called using the same language as Stan, allowing you to use this as a quick reference guide for the look of a given distribution and decide if its an appropriate fit for your data.
Every plot contains the same set of distributions found on the wikipedia page examples of the given distribution.
Documentation on all Stan distributions can be found in Chapter IX of the Stan Reference Guide.
Hope you find this helpful! You can find me on: Github | Twitter | ojhall94-at-gmail-dot-com
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
In the majority (but not all!) of cases, the following holds true:
$\mu$ = the position of the distribution
$\sigma$ = the spread of the distribution.
I have followed the nomenclature used int eh Stan Reference Guide.
$y \sim \textbf{normal}(\textit{mu, sigma})$
def normal(x, mu, sigma):
return (1/np.sqrt(2*np.pi*sigma**2)) * np.exp(-(x - mu)**2/(2*sigma**2))
x = np.linspace(-5,5,100)
plt.plot(x, normal(x, 0, np.sqrt(0.2)))
plt.plot(x, normal(x, 0, np.sqrt(1.0)))
plt.plot(x, normal(x, 0, np.sqrt(5.0)))
plt.plot(x, normal(x, -2, np.sqrt(0.5)))
plt.show()
$y \sim \textbf{exp\_mod\_normal}(\textit{mu, sigma, lambda})$
from scipy.special import erfc
def exp_mod_normal(x, mu, sigma, lambd):
A = (lambd/2) * np.exp((lambd/2)*(2*mu+lambd*sigma**2-2*x))
return A*erfc((mu+lambd*sigma**2-x)/(np.sqrt(2)*sigma))
x = np.linspace(-4,6,100)
plt.plot(x, exp_mod_normal(x, 0, 1., 1.))
plt.plot(x, exp_mod_normal(x, -2, 1., 1.))
plt.plot(x, exp_mod_normal(x, 0, 3., 1.))
plt.plot(x, exp_mod_normal(x, -3, 1., 0.25))
plt.show()
$y \sim \textbf{skew\_normal}(\textit{xi, omega, alpha})$
from scipy.special import erf
def skew_normal(x, xi, omega, alpha):
A = (1/omega*np.sqrt(2*np.pi)) * np.exp(-0.5*((x-xi)/omega)**2)
return A * (1+erf(alpha*((x-xi)/(omega*np.sqrt(2)))))
x = np.linspace(-3,3,100)
plt.plot(x, skew_normal(x, 0., 1., -4.))
plt.plot(x, skew_normal(x, 0., 2., -1.))
plt.plot(x, skew_normal(x, 0., 1., 0.))
plt.plot(x, skew_normal(x, 0., 1., 1.))
plt.plot(x, skew_normal(x, 0., 1., 4.))
plt.show()
$y \sim \textbf{student\_t}(\textit{nu, mu, sigma})$
from scipy.stats import t
def student_t(x, nu, mu, sigma):
return t.pdf(x, nu, mu, sigma)
x = np.linspace(-5,5,100)
plt.plot(x, student_t(x, 1., 0., 1.))
plt.plot(x, student_t(x, 1., 0., 2.))
plt.plot(x, student_t(x, 5., 2., 1.))
plt.plot(x, student_t(x, 1e6, 0., 1.))
plt.show()
$y \sim \textbf{cauchy}(\textit{mu, sigma})$
def cauchy(x, mu, sigma):
return 1/(np.pi*sigma*(1+((x-mu)/sigma)**2))
x = np.linspace(-5,5,100)
plt.plot(x, cauchy(x, 0., 0.5))
plt.plot(x, cauchy(x, 0., 1.))
plt.plot(x, cauchy(x, 0., 2.))
plt.plot(x, cauchy(x, -2., 1.))
plt.show()
$y \sim \textbf{double\_exponential}(\textit{mu, sigma})$
def double_exponential(x, mu, sigma):
return (1./(2.*sigma))*np.exp(-np.abs(x-mu)/sigma)
x = np.linspace(-10,10,100)
plt.plot(x, double_exponential(x, 0.,1))
plt.plot(x, double_exponential(x, 0., 2))
plt.plot(x, double_exponential(x, 0., 4))
plt.plot(x, double_exponential(x, -5., 4))
plt.show()
$y \sim \textbf{logistic}(\textit{mu, sigma})$
def logistic(x, mu, sigma):
return (1./sigma)*np.exp(-(x-mu)/sigma) * (1.+np.exp(-(x-mu)/sigma))**-2
x = np.linspace(-5,20,100)
plt.plot(x, logistic(x, 5.,2.))
plt.plot(x, logistic(x, 9.,3.))
plt.plot(x, logistic(x, 9.,4.))
plt.plot(x, logistic(x, 6.,2.))
plt.plot(x, logistic(x, 2.,1.))
plt.show()
$y \sim \textbf{gumbel}(\textit{mu, beta})$
def gumbel(x, mu, beta):
return (1/beta)*np.exp(-(x-mu)/beta - np.exp(-(x-mu)/beta))
x = np.linspace(-5,20,100)
plt.plot(x, gumbel(x,.5,2.))
plt.plot(x, gumbel(x, 1.,2.))
plt.plot(x, gumbel(x, 1.5,3.))
plt.plot(x, gumbel(x, 3.,4.))
plt.show()
$y \sim \textbf{lognormal}(\textit{mu, sigma})$
def lognormal(x, mu, sigma):
A = (1/np.sqrt(2*np.pi*sigma**2)) * (1/x)
return A * np.exp(-0.5*((np.log(x)-mu)/sigma)**2)
x = np.linspace(0.0001,3,100)
plt.plot(x, lognormal(x,0.,.25))
plt.plot(x, lognormal(x,0.,.5))
plt.plot(x, lognormal(x,0.,1.))
plt.plot(x, lognormal(x,1.,1.))
plt.show()
$y \sim \textbf{chi\_square}(\textit{nu})$
from scipy.stats import chi2
def chi_square(x, nu):
return chi2.pdf(x, nu)
x = np.linspace(0.,8,100)
plt.plot(x, chi_square(x, 1.))
plt.plot(x, chi_square(x, 2.))
plt.plot(x, chi_square(x, 3.))
plt.plot(x, chi_square(x, 6.))
plt.plot(x, chi_square(x, 9.))
plt.show()
$y \sim \textbf{inv\_chi\_square}(\textit{nu})$
from scipy.special import gamma as gammafn
def inv_chi_square(x, nu):
A = (2**(-nu/2)/gammafn(nu/2)) * x**(-nu/2 -1)
return A * np.exp(-1/(2*x))
x = np.linspace(0.0001,1,100)
plt.plot(x, inv_chi_square(x, 1.))
plt.plot(x, inv_chi_square(x, 2.))
plt.plot(x, inv_chi_square(x, 3.))
plt.plot(x, inv_chi_square(x, 6.))
plt.plot(x, inv_chi_square(x, 9.))
plt.show()
$y \sim \textbf{scaled\_inv\_chi\_square}(\textit{nu, sigma})$
from scipy.special import gamma as gammafn
def scaled_inv_chi_square(x, nu, sigma):
A = ((nu/2)**(nu/2) / gammafn(nu/2)) * sigma**nu * x**-(nu/2 +1)
return A * np.exp(-(nu*sigma**2)/(2*x))
x = np.linspace(0.0001,5,100)
plt.plot(x, scaled_inv_chi_square(x, 1., 1.))
plt.plot(x, scaled_inv_chi_square(x, 2., 1.))
plt.plot(x, scaled_inv_chi_square(x, 4., 1.))
plt.plot(x, scaled_inv_chi_square(x, 10., 1.))
plt.show()
$y \sim \textbf{exponential}(\textit{beta})$
def exponential(x, beta):
return beta * np.exp(-beta*x)
x = np.linspace(0.0001,5,100)
plt.plot(x, exponential(x, 0.5))
plt.plot(x, exponential(x, 1.))
plt.plot(x, exponential(x, 1.5))
plt.show()
$y \sim \textbf{gamma}(\textit{alpha, beta})$
from scipy.special import gamma as gammafn
def gamma(x, alpha, beta):
return ((beta**alpha)/gammafn(alpha)) * x**(alpha-1) * np.exp(-beta*x)
x = np.linspace(0.0001, 20, 100)
plt.plot(x, gamma(x, 1., 1./2.))
plt.plot(x, gamma(x, 2., 1./2.))
plt.plot(x, gamma(x, 9., 1./.5))
plt.plot(x, gamma(x, 7.5, 1./1.))
plt.plot(x, gamma(x, .5, 1./1.))
plt.ylim(0.,0.5)
plt.show()
$y \sim \textbf{inv\_gamma}(\textit{alpha, beta})$
from scipy.special import gamma as gammafn
def inv_gamma(x, alpha, beta):
return ((beta**alpha)/gammafn(alpha)) * x**-(alpha+1) * np.exp(-beta*(1/x))
x = np.linspace(0.0001, 3, 100)
plt.plot(x, inv_gamma(x, 1., 1.))
plt.plot(x, inv_gamma(x, 2., 1.))
plt.plot(x, inv_gamma(x, 3., 1.))
plt.plot(x, inv_gamma(x, 3., .5))
plt.show()
$y \sim \textbf{weibull}(\textit{alpha, sigma})$
def weibull(x, alpha, sigma):
return (alpha/sigma)*(x/sigma)**(alpha-1) * np.exp(-(x/sigma)**alpha)
x = np.linspace(0.0001, 2.5, 100)
plt.plot(x, weibull(x, 0.5, 1.))
plt.plot(x, weibull(x, 1., 1.))
plt.plot(x, weibull(x, 1.5, 1.))
plt.plot(x, weibull(x, 5., 1.))
plt.ylim(0., 2.5)
plt.show()
$y \sim \textbf{frechet}(\textit{alpha, sigma})$
def frechet(x, alpha, sigma):
return (alpha/sigma)*(x/sigma)**(-alpha-1) *np.exp(-(x/sigma)**-alpha)
x = np.linspace(0.0001, 4., 100)
plt.plot(x, frechet(x, 1., 1.))
plt.plot(x, frechet(x, 2, 1.))
plt.plot(x, frechet(x, 3, 1.))
plt.plot(x, frechet(x, 1, 2.))
plt.plot(x, frechet(x, 2, 2.))
plt.plot(x, frechet(x, 3, 2.))
plt.show()
$y \sim \textbf{rayleigh}(\textit{sigma})$
def rayleigh(x, sigma):
return (x/sigma**2)*np.exp(-x**2/(2*sigma**2))
x = np.linspace(0., 10., 100)
plt.plot(x, rayleigh(x, .5))
plt.plot(x, rayleigh(x, 1.))
plt.plot(x, rayleigh(x, 2.))
plt.plot(x, rayleigh(x, 3.))
plt.plot(x, rayleigh(x, 4.))
plt.show()
$y \sim \textbf{pareto}(\textit{y_min, alpha})$
def pareto(x, x_min, alpha):
return (alpha*x_min**alpha)/(x**(alpha+1))
x = np.linspace(1., 5., 100)
plt.plot(x, pareto(x, 1., 1.))
plt.plot(x, pareto(x, 1., 2.))
plt.plot(x, pareto(x, 1., 3.))
plt.plot(x, pareto(x, 3., 2.))
plt.ylim(0.,3.)
plt.show()
Note: The low end of the distribution ends at the given value of y_min.
$y \sim \textbf{pareto\_type\_2}(\textit{mu, lambda, alpha})$
def pareto_type_2(x, mu, lambd, alpha):
return (alpha/lambd)*(1+((x-mu)/lambd))**(-(alpha+1))
x = np.linspace(0., 6., 100)
plt.plot(x, pareto_type_2(x, 0., 1., 2.))
plt.plot(x, pareto_type_2(x, 0., 2., 2.))
plt.plot(x, pareto_type_2(x, 0., 4., 1.))
plt.plot(x, pareto_type_2(x, 0., 6., 1.))
plt.show()
$y \sim \textbf{beta}(\textit{alpha, beta})$
from scipy.special import beta as betafn
def beta(theta, alpha, beta):
return (1/betafn(alpha, beta))*x**(alpha-1) * (1-x)**(beta-1)
x = np.linspace(0.0001, 0.9999, 100)
plt.plot(x, beta(x, .5, .5))
plt.plot(x, beta(x, 5., 1.))
plt.plot(x, beta(x, 1., 3.))
plt.plot(x, beta(x, 2., 2.))
plt.plot(x, beta(x, 2., 5.))
plt.ylim(0., 2.5)
plt.show()
Everything else in the reference guide is slightly above my ability and need to understand right now. Feel free to add them in though, if you find yourself using them!