Data

rng = np.random.default_rng(seed = 1234)
cl1 = rng.multivariate_normal([-2,-2], [[1,-0.5],[-0.5,1]], size=100)
cl2 = rng.multivariate_normal([1,0], [[1,0],[0,1]], size=150)
cl3 = rng.multivariate_normal([3,2], [[1,-0.7],[-0.7,1]], size=200)
pts = np.concatenate((cl1,cl2,cl3))

k-means clustering

from scipy.cluster.vq import kmeans
ctr, dist = kmeans(pts, 3)

ctr

array([[ 2.86399,  1.95401],
       [-2.03957, -1.85662],
       [ 0.91124, -0.18724]])

dist

np.float64(1.2209235923868729)

cl1.mean(axis=0)

array([-2.00475, -1.87276])

cl2.mean(axis=0)

array([1.03849, 0.01417])

cl3.mean(axis=0)

array([2.94642, 2.02514])

k-means distortion plot

The mean (non-squared) Euclidean distance between the observations passed and the centroids generated.

ks = range(1,8)
dists = [kmeans(pts, k)[1] for k in ks]

np.array(dists).reshape(-1)

array([2.54703, 1.57178, 1.22092, 1.04749, 0.95422, 0.88013, 0.8176 ])

Assigning new points with vq()

Once we have centroids from kmeans(), we can assign new points to clusters using vq() (vector quantization).

from scipy.cluster.vq import vq

new_pts = np.array([[0, 0], [-3, -3], [4, 3]])
labels, dists = vq(new_pts, ctr)

labels

array([2, 1, 0], dtype=int32)

dists

array([0.93028, 1.49323, 1.54422])

Basic functions

For general numeric integration in 1D we use scipy.integrate.quad(), which takes as arguments: the function to be integrated and then the lower and upper bounds of the integral.

from scipy.integrate import quad

quad(lambda x: x, 0, 1)

(0.5, 5.551115123125783e-15)

quad(np.sin, 0, np.pi)

(2.0, 2.220446049250313e-14)

quad(np.sin, 0, 2*np.pi)

(2.221501482512777e-16, 4.3998892617846e-14)

quad(np.exp, 0, 1)

(1.7182818284590453, 1.9076760487502457e-14)

Normal PDF

The PDF for a normal distribution is given by,

\[ f(x) = \frac{1}{\sigma \sqrt{2 \pi}} \exp\left(-\frac{1}{2} \left(\frac{x-\mu}{\sigma}\right)^2 \right) \]

def norm_pdf(x, μ, σ):
  return (1/(σ * np.sqrt(2*np.pi))) * np.exp(-0.5 * ((x - μ)/σ)**2)

norm_pdf(0, 0, 1)

np.float64(0.3989422804014327)

norm_pdf(np.inf, 0, 1)

np.float64(0.0)

norm_pdf(-np.inf, 0, 1)

np.float64(0.0)

Checking the PDF

We can check that we’ve implemented a valid pdf by integrating from \(-\infty\) to \(\infty\),

quad(norm_pdf, -np.inf, np.inf)

TypeError: norm_pdf() missing 2 required positional arguments: 'μ' and 'σ'

quad(lambda x: norm_pdf(x, 0, 1), -np.inf, np.inf)

(0.9999999999999997, 1.0178191380347127e-08)

quad(lambda x: norm_pdf(x, 17, 12), -np.inf, np.inf)

(1.0000000000000002, 4.113136862574909e-09)

Truncated normal PDF

\[ f(x) = \begin{cases} \frac{c}{\sigma \sqrt{2 \pi}} \exp\left(-\frac{1}{2} \left(\frac{x-\mu}{\sigma}\right)^2 \right), & \text{for } a \leq x \leq b \\ 0, & \text{otherwise.} \\ \end{cases} \]

def trunc_norm_pdf(x, μ=0, σ=1, a=-np.inf, b=np.inf):
  if (b < a):
      raise ValueError("b must be greater than a")

  x = np.asarray(x)
  scalar_input = x.ndim == 0
  x = np.atleast_1d(x)

  full_pdf = (1/(σ * np.sqrt(2*np.pi))) * np.exp(-0.5 * ((x - μ)/σ)**2)
  full_pdf[(x < a) | (x > b)] = 0

  return full_pdf[0] if scalar_input else full_pdf

Testing trunc_norm_pdf

trunc_norm_pdf(0, a=-1, b=1)

np.float64(0.3989422804014327)

trunc_norm_pdf(2, a=-1, b=1)

np.float64(0.0)

trunc_norm_pdf(-2, a=-1, b=1)

np.float64(0.0)

trunc_norm_pdf([-2,1,0,1,2], a=-1, b=1)

array([0.     , 0.24197, 0.39894, 0.24197, 0.     ])

quad(lambda x: trunc_norm_pdf(x, a=-1, b=1), -np.inf, np.inf)

(0.682689492137086, 2.0147661317082566e-11)

quad(lambda x: trunc_norm_pdf(x, a=-3, b=3), -np.inf, np.inf)

(0.9973002039367396, 7.451935936375609e-09)

Fixing trunc_norm_pdf

def trunc_norm_pdf(x, μ=0, σ=1, a=-np.inf, b=np.inf):
  if (b < a):
      raise ValueError("b must be greater than a")

  x = np.asarray(x)
  scalar_input = x.ndim == 0
  x = np.atleast_1d(x)

  nc = 1 / quad(lambda x: norm_pdf(x, μ, σ), a, b)[0]

  full_pdf = nc * (1/(σ * np.sqrt(2*np.pi))) * np.exp(-0.5 * ((x - μ)/σ)**2)
  full_pdf[(x < a) | (x > b)] = 0

  return full_pdf[0] if scalar_input else full_pdf

Multivariate normal

\[ f(\bf{x}) = \det{(2\pi\Sigma)}^{-1/2} \exp{\left(-\frac{1}{2} (\bf{x}-\mu)^T \Sigma^{-1}(\bf{x}-\mu) \right)} \]

def mv_norm(x, μ, Σ):
  x = np.asarray(x)
  μ = np.asarray(μ)
  Σ = np.asarray(Σ)
  
  return ( np.linalg.det(2*np.pi*Σ)**(-0.5) * 
           np.exp(-0.5 * (x - μ).T @ np.linalg.solve(Σ, (x-μ)) ) )

norm_pdf(0,0,1)

np.float64(0.3989422804014327)

mv_norm([0], [0], [[1]])

np.float64(0.3989422804014327)

mv_norm([0,0], [0,0], [[1,0],[0,1]])

np.float64(0.15915494309189535)

mv_norm([0,0,0], [0,0,0], 
        [[1,0,0],[0,1,0],[0,0,1]])

np.float64(0.06349363593424098)

2d & 3d numerical integration

These are supported by dblquad() and tplquad() respectively (see nquad() for higher dimensions).

from scipy.integrate import dblquad, tplquad

dblquad(lambda y, x: mv_norm([x,y], [0,0], np.identity(2)),
        a=-np.inf, b=np.inf,
        gfun=lambda x: -np.inf,   hfun=lambda x: np.inf)

(1.0000000000000322, 1.315012783659768e-08)

tplquad(lambda z, y, x: mv_norm([x,y,z], [0,0,0], np.identity(3)),
        a=0, b=np.inf,
        gfun=lambda x:   0, hfun=lambda x:   np.inf,
        qfun=lambda x,y: 0, rfun=lambda x,y: np.inf)

(0.12500000000036066, 1.4697203688869196e-08)

Scalar function minimization

def f(x):
    return x**4 + 3*(x-2)**3 - 15*(x)**2 + 1

from scipy.optimize import minimize_scalar
minimize_scalar(f)

 message: 
          Optimization terminated successfully;
          The returned value satisfies the termination criteria
          (using xtol = 1.48e-08 )
 success: True
     fun: -803.3955308825884
       x: -5.528801125219663
     nit: 11
    nfev: 16

Results

res = minimize_scalar(f)

type(res)

scipy.optimize._optimize.OptimizeResult

dir(res)

['fun', 'message', 'nfev', 'nit', 'success', 'x']

res.fun

np.float64(-803.3955308825884)

print(res.message)

Optimization terminated successfully;
The returned value satisfies the termination criteria
(using xtol = 1.48e-08 )

res.nfev

16

res.nit

11

res.success

True

res.x

np.float64(-5.528801125219663)

More details

from scipy.optimize import show_options
show_options(solver="minimize_scalar")

brent
=====

Options
-------
maxiter : int
    Maximum number of iterations to perform.
xtol : float
    Relative error in solution `xopt` acceptable for convergence.
disp : int, optional
    If non-zero, print messages.

    ``0`` : no message printing.

    ``1`` : non-convergence notification messages only.

    ``2`` : print a message on convergence too.

    ``3`` : print iteration results.

Notes
-----
Uses inverse parabolic interpolation when possible to speed up
convergence of golden section method.

bounded
=======

Options
-------
maxiter : int
    Maximum number of iterations to perform.
disp: int, optional
    If non-zero, print messages.

    ``0`` : no message printing.

    ``1`` : non-convergence notification messages only.

    ``2`` : print a message on convergence too.

    ``3`` : print iteration results.

xatol : float
    Absolute error in solution `xopt` acceptable for convergence.

golden
======

Options
-------
xtol : float
    Relative error in solution `xopt` acceptable for convergence.
maxiter : int
    Maximum number of iterations to perform.
disp: int, optional
    If non-zero, print messages.

    ``0`` : no message printing.

    ``1`` : non-convergence notification messages only.

    ``2`` : print a message on convergence too.

    ``3`` : print iteration results.

Local minima

def f(x):
  return -np.sinc(x-5)

res = minimize_scalar(f); res

 message: 
          Optimization terminated successfully;
          The returned value satisfies the termination criteria
          (using xtol = 1.48e-08 )
 success: True
     fun: -0.049029624014074166
       x: -1.4843871263953001
     nit: 10
    nfev: 14

Random starts

rng = np.random.default_rng(seed=1234)

lower = rng.uniform(-20, 20, 100)
upper = lower + 1

sols = [minimize_scalar(f, bracket=(l,u)) 
        for l,u in zip(lower, upper)]
funs = [sol.fun for sol in sols]

best = sols[np.argmin(funs)]
best

 message: 
          Optimization terminated successfully;
          The returned value satisfies the termination criteria
          (using xtol = 1.48e-08 )
 success: True
     fun: -1.0
       x: 5.0000000006185585
     nit: 8
    nfev: 11

Distributions

Implements classes for most continuous and discrete distributions,

• rvs - Random Variates

• pdf - Probability Density Function

• cdf - Cumulative Distribution Function

• sf - Survival Function (1-CDF)

• ppf - Percent Point Function (Inverse of CDF)

• isf - Inverse Survival Function (Inverse of SF)

• stats - Return mean, variance, (Fisher’s) skew, or (Fisher’s) kurtosis

• moment - non-central moments of the distribution

Basic usage

from scipy.stats import norm, gamma, binom, uniform

norm().rvs(size=5)

array([-0.27378, -0.23604, -0.02471,  0.68152,  0.12815])

uniform.pdf([0,0.5,1,2])

array([1., 1., 1., 0.])

binom.mean(n=10, p=0.25)

np.float64(2.5)

binom.median(n=10, p=0.25)

np.float64(2.0)

norm().stats()

(np.float64(0.0), np.float64(1.0))

gamma(a=1,scale=1).stats(moments="mvsk")

(np.float64(1.0), np.float64(1.0), np.float64(2.0), np.float64(6.0))

Freezing

Model parameters can be passed to any of the methods directly, or a distribution can be constructed using a specific set of parameters, which is known as freezing.

norm_rv = norm(loc=-1, scale=3)
norm_rv.median()

np.float64(-1.0)

unif_rv = uniform(loc=-1, scale=2)
unif_rv.cdf([-2,-1,0,1,2])

array([0. , 0. , 0.5, 1. , 1. ])

unif_rv.rvs(5)

array([-0.58057,  0.46872, -0.8764 ,  0.83672, -0.16421])

g = gamma(a=2, loc=0, scale=1.2)

plt.figure(figsize=(3, 2))
x = np.linspace(0, 10, 100)
plt.plot(x, g.pdf(x), "k-")
plt.axvline(x=g.mean(), c="r")
plt.axvline(x=g.median(), c="b")

MLE

Maximum likelihood estimation is possible via the fit() method,

x = norm.rvs(loc=2.5, scale=2, size=1000, random_state=1234)

norm.fit(x)

(np.float64(2.5314811643075235), np.float64(1.946132398754459))

norm.fit(x, loc=2.5) # provide a guess for the parameter

(np.float64(2.5314811643075235), np.float64(1.946132398754459))

x = gamma.rvs(a=2.5, size=1000)
gamma.fit(x) # shape, loc, scale

(np.float64(2.589173185346138),
 np.float64(-0.00011290008283903167),
 np.float64(0.9676759000682381))

y = gamma.rvs(a=2.5, loc=-1, scale=2, size=1000)
gamma.fit(y) # shape, loc, scale

(np.float64(2.7823124612208137),
 np.float64(-1.1090393239477736),
 np.float64(1.8785702919274456))

scipy.special

Provides implementations of many special mathematical functions commonly used in statistics, physics, and engineering.

Gamma and Beta functions

from scipy.special import gamma, gammaln, beta, betaln

gamma(5)  # (n-1)! for integers

np.float64(24.0)

gamma(0.5)  # sqrt(pi)

np.float64(1.7724538509055159)

np.sqrt(np.pi)

np.float64(1.7724538509055159)

gammaln(100)  # log(gamma(100))

np.float64(359.1342053695754)

beta(2, 3)  # B(a,b) = gamma(a)*gamma(b)/gamma(a+b)

np.float64(0.08333333333333333)

gamma(2) * gamma(3) / gamma(5)

np.float64(0.08333333333333333)

Combinatorial functions

from scipy.special import factorial, comb, perm

factorial(5)

np.float64(120.0)

factorial(np.arange(6))

array([  1.,   1.,   2.,   6.,  24., 120.])

comb(10, 3)  # 10 choose 3

np.float64(120.0)

perm(10, 3)  # 10 permute 3

np.float64(720.0)

Error function

The error function \(\text{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt\),

from scipy.special import erf, erfc

erf(0)

np.float64(0.0)

erf(np.inf)

np.float64(1.0)

erfc(0)  # 1 - erf(x)

np.float64(1.0)

Logistic functions

from scipy.special import expit, logit, softmax

expit(0)  # 1 / (1 + exp(-x))

np.float64(0.5)

expit([-np.inf, 0, np.inf])

array([0. , 0.5, 1. ])

logit(0.5)  # log(p / (1-p))

np.float64(0.0)

logit([0.1, 0.5, 0.9])

array([-2.19722,  0.     ,  2.19722])

softmax([1, 2, 3])

array([0.09003, 0.24473, 0.66524])

softmax([1, 2, 3]).sum()

np.float64(0.9999999999999999)

softmax([0, 0, 0])

array([0.33333, 0.33333, 0.33333])

NumPy vs SciPy linalg

Both numpy.linalg and scipy.linalg provide linear algebra routines, but SciPy’s version is more comprehensive.

Basic operations comparison

import scipy.linalg

A = np.array([[1, 2], [3, 4]])
b = np.array([5, 6])

np.linalg.solve(A, b)

array([-4. ,  4.5])

np.linalg.inv(A)

array([[-2. ,  1. ],
       [ 1.5, -0.5]])

scipy.linalg.solve(A, b)

array([-4. ,  4.5])

scipy.linalg.inv(A)

array([[-2. ,  1. ],
       [ 1.5, -0.5]])

SciPy-only features

P, L, U = scipy.linalg.lu(A)

L

array([[1.     , 0.     ],
       [0.33333, 1.     ]])

U

array([[3.     , 4.     ],
       [0.     , 0.66667]])

scipy.linalg.expm(A)

array([[ 51.96896,  74.73656],
       [112.10485, 164.0738 ]])

scipy.linalg.logm(
  scipy.linalg.expm(A)
)

array([[1., 2.],
       [3., 4.]])

Dense vs Sparse matrices

For matrices with many zero entries, sparse representations are more memory-efficient and can be much faster.

from scipy import sparse

dense = np.array([
    [1, 0, 0, 0],
    [0, 2, 0, 0],
    [0, 0, 3, 0],
    [0, 0, 0, 4]
])
dense.nbytes

128

sp = sparse.csr_matrix(dense)
sp

<Compressed Sparse Row sparse matrix of dtype 'int64'
    with 4 stored elements and shape (4, 4)>

sp.data

array([1, 2, 3, 4])

Sparse matrix formats

Creating sparse matrices

row = [0, 1, 2, 2]
col = [0, 1, 0, 2]
data = [1, 2, 3, 4]
sp = sparse.coo_matrix(
  (data, (row, col)), shape=(3, 3)
)
sp.toarray()

array([[1, 0, 0],
       [0, 2, 0],
       [3, 0, 4]])

sparse.eye(4, format="csr").toarray()

array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]])

sparse.random(3, 3, density=0.3).toarray()

array([[0.     , 0.     , 0.     ],
       [0.     , 0.     , 0.22213],
       [0.     , 0.30699, 0.68624]])

sparse.diags([1., 2., 3., 4.]).toarray()

array([[1., 0., 0., 0.],
       [0., 2., 0., 0.],
       [0., 0., 3., 0.],
       [0., 0., 0., 4.]])

Sparse linear algebra

scipy.sparse.linalg provides solvers optimized for sparse matrices.

from scipy.sparse.linalg import spsolve, eigs

n = 1000
diagonals = [np.ones(n-1), -2*np.ones(n), np.ones(n-1)]
A_sparse = sparse.diags(diagonals, [-1, 0, 1], format="csr")
b = np.ones(n)

x = spsolve(A_sparse, b)
x[:5]

array([ -500.,  -999., -1497., -1994., -2490.])

vals, vecs = eigs(A_sparse, k=3)
vals

array([-3.99999+0.j, -3.99996+0.j, -3.99991+0.j])

When to use sparse?

General rule - use sparse when density < 10% and matrix is large.