stan ingredients

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

In [1]:
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.

Unbounded Continuous Distributions

Normal Distribution

$y \sim \textbf{normal}(\textit{mu, sigma})$

In [2]:
def normal(x, mu, sigma):
    return (1/np.sqrt(2*np.pi*sigma**2)) * np.exp(-(x - mu)**2/(2*sigma**2))
In [3]:
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()

Exponentially Modified Normal Distribution

$y \sim \textbf{exp\_mod\_normal}(\textit{mu, sigma, lambda})$

In [4]:
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))
In [5]:
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()

Skew Normal Distribution

$y \sim \textbf{skew\_normal}(\textit{xi, omega, alpha})$

In [6]:
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)))))
In [7]:
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()

Student-t Distribution

$y \sim \textbf{student\_t}(\textit{nu, mu, sigma})$

In [8]:
from scipy.stats import t

def student_t(x, nu, mu, sigma):
    return t.pdf(x, nu, mu, sigma)
In [9]:
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()

Cauchy (Lorentzian) Distribution

$y \sim \textbf{cauchy}(\textit{mu, sigma})$

In [10]:
def cauchy(x, mu, sigma):
    return 1/(np.pi*sigma*(1+((x-mu)/sigma)**2))
In [11]:
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()

Double Exponential (Laplace) Distribution

$y \sim \textbf{double\_exponential}(\textit{mu, sigma})$

In [12]:
def double_exponential(x, mu, sigma):
    return (1./(2.*sigma))*np.exp(-np.abs(x-mu)/sigma)
In [13]:
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()

Logistic Distribution

$y \sim \textbf{logistic}(\textit{mu, sigma})$

In [14]:
def logistic(x, mu, sigma):
    return (1./sigma)*np.exp(-(x-mu)/sigma) * (1.+np.exp(-(x-mu)/sigma))**-2
In [15]:
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()

Gumbel Distribution

$y \sim \textbf{gumbel}(\textit{mu, beta})$

In [16]:
def gumbel(x, mu, beta):
    return (1/beta)*np.exp(-(x-mu)/beta - np.exp(-(x-mu)/beta))
In [17]:
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()

Positive Continuous Distributions

Lognormal Distribution

$y \sim \textbf{lognormal}(\textit{mu, sigma})$

In [18]:
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)
In [19]:
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()

Chi-Square Distribution

$y \sim \textbf{chi\_square}(\textit{nu})$

In [20]:
from scipy.stats import chi2
def chi_square(x, nu):
    return chi2.pdf(x, nu)
In [21]:
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()

Inverse Chi-Square Distribution

$y \sim \textbf{inv\_chi\_square}(\textit{nu})$

In [22]:
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))
In [23]:
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()

Scaled Inverse Chi-Square Distribution

$y \sim \textbf{scaled\_inv\_chi\_square}(\textit{nu, sigma})$

In [24]:
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))
In [25]:
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()

Exponential Distribution

$y \sim \textbf{exponential}(\textit{beta})$

In [26]:
def exponential(x, beta):
    return beta * np.exp(-beta*x)
In [27]:
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()

Gamma Distribution

$y \sim \textbf{gamma}(\textit{alpha, beta})$

In [28]:
from scipy.special import gamma as gammafn
def gamma(x, alpha, beta):
    return ((beta**alpha)/gammafn(alpha)) * x**(alpha-1) * np.exp(-beta*x)
In [29]:
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()

Inverse Gamma Distribution

$y \sim \textbf{inv\_gamma}(\textit{alpha, beta})$

In [30]:
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))
In [31]:
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()

Weibull Distribution

$y \sim \textbf{weibull}(\textit{alpha, sigma})$

In [32]:
def weibull(x, alpha, sigma):
    return (alpha/sigma)*(x/sigma)**(alpha-1) * np.exp(-(x/sigma)**alpha)
In [33]:
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()

Frechet Distribution

$y \sim \textbf{frechet}(\textit{alpha, sigma})$

In [34]:
def frechet(x, alpha, sigma):
    return (alpha/sigma)*(x/sigma)**(-alpha-1) *np.exp(-(x/sigma)**-alpha)
In [35]:
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()

Non-negative Continuous Distributions

Rayleigh Distribution

$y \sim \textbf{rayleigh}(\textit{sigma})$

In [36]:
def rayleigh(x, sigma):
    return (x/sigma**2)*np.exp(-x**2/(2*sigma**2))
In [37]:
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()

Positive Lower Bounded Probabilities

Pareto Distribution

$y \sim \textbf{pareto}(\textit{y_min, alpha})$

In [38]:
def pareto(x, x_min, alpha):
    return (alpha*x_min**alpha)/(x**(alpha+1))
In [39]:
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.

Pareto Type 2 (Lomax) Distribution

$y \sim \textbf{pareto\_type\_2}(\textit{mu, lambda, alpha})$

In [40]:
def pareto_type_2(x, mu, lambd, alpha):
    return (alpha/lambd)*(1+((x-mu)/lambd))**(-(alpha+1))
In [41]:
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()

Continuous Distributions on [0, 1]

Beta Distribution

$y \sim \textbf{beta}(\textit{alpha, beta})$

In [42]:
from scipy.special import beta as betafn
def beta(theta, alpha, beta):
    return (1/betafn(alpha, beta))*x**(alpha-1) * (1-x)**(beta-1)
In [43]:
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()

Thats all folks!

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!