Distributions

The bernoulli distribution represents the frequency of successes in a given number of trials. It is just a special case of a binomial distribution (the expression of which we opened this paper) where the number of observations \(N\) is \(1\) and, thus, the possible outcomes are \(y \in [0,1]\). Let’s simulate coin tosses in R by randomly sampling from possible events in the distribution space,

set.seed(TRUE)  # for reproducibility

flips <- sample(x = c("HH", "TT", "HT", "TH"),    # sample space
                size = 1000,                      # number of trials
                prob = c(0.25, 0.25, 0.25, 0.25), # probability of each event
                replace = TRUE)                   # replace the coins each time

flips <- factor(x = flips, labels = c("HH", "TT", "HT", "TH"))
summary(flips)

##  HH  TT  HT  TH 
## 248 276 232 244

Using our presumed probability of obtaining two heads on any given flip as .25, let’s calculate the probability of two heads occuring exactly 248 times out of our 1000 draws, first expressly coding the binomial function, and then using R’s built-in function.

p1 <- choose(1000, sum(flips == 'HH')) * .25 ^ 248 * (1 - .25) ^ (1000 - sum(flips == 'HH'))
p1

## [1] 0.02889194

p2 <- dbinom(sum(flips=='HH'), 1000, .25)
p2

## [1] 0.02889194

We can test to see that both methods are equal:

all.equal(p1, p2)

## [1] TRUE

R offers several functions for common distributions: a probability mass function (PMF) for discrete variables or probability density function (PDF) for continuous variables, a continuous distribution function (CDF), and a function to draw random values from the distribution. R’s naming convention for these functions begin with a letter describing the function type — p for a CDF, d for a PMF/PDF, and r for a random generation function — followed by letters identifying the distribution. Thus, the above dbinom() provides the PMF for the binomial. Let’s apply these and review them in a graph. For a binomial distribution with a 50 percent probability of success in each of 40 observations of that many draws, for example, we can graph its PMF, CDF, and a histogram of these random draws,

require(ggplot2)                               # load plotting library

x <- seq(1, 40, 1)                             # create a sequence of values

df <- data.frame(value = x,                    # organize values, cdf, pdf
                 cdf = pbinom(x, length(x), .5),
                 pmf = dbinom(x, length(x), .5),
                 draws = rbinom(length(x), length(x), .5))

ggplot(data = df) +                            # create plot object

  # plot density of draws
  geom_histogram(aes(x = draws,                
                     y = ..density..), 
                 position = 'identity',          
                 binwidth = 1,
                 color = 'white',
                 fill = 'gray70') +
  
  # plot the CDF
  geom_line(aes(x = value,                     
                y = cdf),
            color = '#91A5AF') + 
  
  # plot the PMF
  geom_line(aes(x = value,                     
                y = pmf)) +
  
  # remove default grid lines
  theme(panel.grid.major = element_blank(),    
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) +
  
  # add labels to plot elements
  annotate('text', x = 25, y = 0.8, label = 'CDF', color = '#91A5AF') +
  annotate('text', x = 30, y = 0.1, label = 'PMF', color = 'black') +
  annotate('text', x = 21, y = 0.25, label = 'DRAWS', color = 'gray70') +
  
  # add labels to axes
  labs(x = "Successes of 40\ndraws in 40 trials",                            
       y = "Probability")

(Recall that a PDF, unlike a PMF, does not represent probability; it represents the slope of the CDF.) Other common (and interrelated) distributions encountered include the multinomial dmultinom(), normal or gaussian dnorm(), cauchy dcauchy(), student-t dt(), poisson dpois(), negative binomial dnbinom(), exponential dexp(), gamma dgamma(), beta dbeta(). The LKJ correlation matrix distribution is dlkjcorr().

Some probability is reviewed in chapter three of Deep Learing. More briefly, we get a quick reference on probability from Chen, et al. Probability Cheatsheet, (2015), which closely follows Blitzstein’s Introduction to Probability (CRC Press 2014). A clear and concise statement on probability is found in, believe it or not, Feynman’s Chapter 6 Probability in Lectures on Physics, Volume I Millennium Edition (Basic Books 2010). Another helpful review is found in Murphy’s Machine Learning, A Probablistic Perspective, Chapter 2 Probability (MIT Press 2012).

Regression

Borrowing from Gelman and Hill’s Data Analysis Using Regression and Multilevel/Heirarchical Models (Cambridge University Press 2007), linear “regression can be used to represent relationships between variables,” and “is a method that summarizes how the average values of a numerical outcome variable vary over subpopulations defined by linear functions of predictors.” With mathematical notation, we can represent Gelman’s description as \(y_i = \boldsymbol{X}_i \boldsymbol{\beta}+\epsilon = \beta_1\boldsymbol{X}_{i1}+\dots+\beta_k\boldsymbol{X}_{ik}+\epsilon\), for \(i=1, \dots, n\), where errors \(\epsilon_i\) are normally distributed with mean \(0\) and standard deviation \(\sigma\).

# prepare data and linear model
# -----------------------------------------------------------------------------
# save toy baseball data into object
df <- subset(plyr::baseball, year == 2007 & ab > 100)

# use basic linear model
mod <- lm(hr ~ ab, data = df)

# add fitted values to data
df <- transform(df, Fitted = fitted(mod))

# Create plot for SSR
# -----------------------------------------------------------------------------
SSR <- ggplot(df, aes(x = ab, y = hr)) + 

  # plot yhat - ybar
  geom_segment(aes(x = ab, 
                   y = mean(hr),
                   xend = ab, 
                   yend = Fitted),
               size = .2,
               color = '#888888') +

  # plot observations as points
  geom_point(color = 'gray85') +

  # plot linear regression line
  geom_smooth(method='lm', 
              formula = y ~ x,
              se = FALSE,
              color = 'black',
              size = 0.5) + 
  
  # include origin for proper reference point in plot and to show scale
  scale_x_continuous(limits = c(80, 550)) + 
  scale_y_continuous(limits = c(0, 40)) +
  
  # plot mean of samples ybar
  geom_hline(aes(yintercept = mean(hr)),
             size = .5,
             color = 'black') +
  
  # remove all default decoration                
  theme_void() +
  
  # add annotations for title, yhat, ybar
  annotate("text", x = 100, y = 30,
           label = 'italic(hat(y))[i]-italic(bar(y))',
           parse = T, hjust = 0) +
  
  annotate("text", x = 500, y= 22, 
           label = 'bolditalic(hat(y))', 
           parse = T, hjust = 0) +
  
  annotate("text", x = 100, y= 14, 
           label = 'italic(bar(y))', 
           parse = T, hjust = 0)

# Create plot for SST
# -----------------------------------------------------------------------------
SST <- ggplot(df, aes(x = ab, y = hr)) + 

  # plot y - ybar
  geom_segment(aes(x = ab, 
                   y = hr,
                   xend = ab, 
                   yend = mean(hr)),
               size = .2,
               color = '#888888') +

  # plot observations as points
  geom_point(color = 'gray30') +

  # plot linear regression line
  geom_smooth(method='lm', 
              formula = y ~ x,
              se = FALSE,
              color = 'gray85',
              size = 0.5) + 
  
  # include origin for proper reference point in plot and to show scale
  scale_x_continuous(limits = c(80, 550)) + 
  scale_y_continuous(limits = c(0, 40)) +
  
  # plot mean of samples ybar
  geom_hline(aes(yintercept = mean(hr)),
             size = .5,
             color = 'black') +
  
  # remove all default decoration                
  theme_void() +
  
  # add annotations for title, yhat, ybar
  annotate("text", x = 100, y = 30,
           label = 'italic(y)[i]-italic(bar(y))',
           parse = T, hjust = 0) +

  annotate("text", x = 500, y= 22, 
           label = 'bolditalic(hat(y))', 
           parse = T, hjust = 0) +
  
  annotate("text", x = 100, y= 14, 
           label = 'italic(bar(y))', 
           parse = T, hjust = 0)

# create plot for epsilon
# -----------------------------------------------------------------------------
epsilon <- ggplot(df, aes(x = ab, y = hr)) + 

  # draw y - yhat
  geom_segment(aes(x = ab, 
                   y = hr,
                   xend = ab, 
                   yend = Fitted),
               size = .2,
               color = '#888888') +

  # plot observations as points
  geom_point(color = 'gray30') +

  # plot linear regression line
  geom_smooth(method='lm', 
              formula = y ~ x,
              se = FALSE,
              color = 'black',
              size = 0.5) + 
  
  # include origin for proper reference point in plot and to show scale
  scale_x_continuous(limits = c(80, 550)) + 
  scale_y_continuous(limits = c(0, 40)) +
  
  # plot mean of samples ybar
  geom_hline(aes(yintercept = mean(hr)),
             size = .3,
             color = 'gray85') +
                  
  # remove all default decoration                
  theme_void() +
  
  # add annotations for title, yhat, ybar
  annotate("text", x = 100, y = 30,
           label = 'italic(y)[i]-italic(hat(y))[i]',
           parse = T, hjust = 0) +

  annotate("text", x = 500, y= 22, 
           label = 'bolditalic(hat(y))', 
           parse = T, hjust = 0) +
  
  annotate("text", x = 100, y= 14, 
           label = 'italic(bar(y))', 
           parse = T, hjust = 0)

# plot SSR, SST, and epsilon together
# -----------------------------------------------------------------------------
require(grid)
require(gridExtra)

grid.arrange(SSR, SST, epsilon, ncol = 3)

Gradient descent

Let’s estimate our linear regression using gradient descent for the gradient \(\nabla f(\theta)\). A gradient as the analogue of the first derivative for functions of vectors. The gradient of a function \(f\) is the matrix of partial derivatives, taken generally as one or more partial derivatives,

\[ \nabla_\theta f(\theta) \in \mathbb{R}^{m \times n} = \left[\begin{matrix} \frac{\partial f(\theta)}{\partial \theta_{11}} & \dots & \frac{\partial f(\theta)}{\partial \theta_{1n}} \\ \vdots & \ddots & \vdots \\ \frac{\partial f(\theta)}{\partial \theta_{m1}} & \dots & \frac{\partial f(\theta)}{\partial \theta_{mn}} \end{matrix}\right] \] and as specific to linear regression,

\[ \nabla f(\theta) = \frac{1}{N}(\boldsymbol{y}^T - \theta \boldsymbol{X}^T)\boldsymbol{X} \] Our iterative procedure is simply,

\[ \text{Repeat until convergence:}\\ \theta := \theta - \alpha \sum{(\boldsymbol{y}^T - \theta \boldsymbol{X}^T)\boldsymbol{X}} \] Let’s code an example in R,

# generate matrix with intercept and random variable
X <- cbind(1, matrix(runif(1000,-5, 5)))
y <- X[,2] + rnorm(1000) + 3

# alpha and threshold manually chosen for toy data
gradient_descent <- function(X, y, alpha = 0.01, threshold = .0001) {
  
  # initialize variables
  gradient <- Inf
  i <- 1
  
  # initialize coefficients
  theta <- matrix(rep(0, dim(X)[2]), nrow = 2)

  # intermediate values
  cost_history <- double()
  theta_history <- list()

  # Repeat until convergence
  while(abs(max(gradient)) > threshold) {
    
    # update theta with gradient each iteration
    gradient <- ( t(X) %*% (X %*% theta - y) / length(y) )
    theta <- theta - alpha * gradient
  
    # store theta and cost (sse) each iteration
    cost_history[i] <- sum((X %*% theta - y) ^ 2) / (2 * length(y))
    theta_history[[i]] <- theta
    i <- i + 1
  }
  
  cost_history <- as.array(cost_history)
  theta_history <- matrix(unlist(theta_history), ncol = 2, byrow = T)
  return(list(cost = cost_history,
              theta = theta_history))
}

gd <- gradient_descent(X, y)

Plotting each iteration of \(\theta\) from our algorithm, along with the data,

ggplot() + 
  
  # plot simulated data
  geom_point(aes(x = X[,2], y = y), color = 'gray90') +

  # plot intermediate gradients
  geom_abline(aes(slope = gd$theta[,2], 
                  intercept = gd$theta[,1]),
              color = 'gray', 
              alpha = .3) +

  # plot starting gradient
  geom_abline(aes(slope = gd$theta[1, 2], 
                  intercept = gd$theta[1, 1]),
              color = 'gray50') +

  # plot final gradient
  geom_abline(aes(slope = gd$theta[nrow(gd$theta),2], 
                  intercept = gd$theta[nrow(gd$theta),1]),
              color = '#000000') +
  
  # clean plot
  theme(panel.grid.major = element_blank(),    
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_blank(),
        plot.title = element_text(hjust = 0.5)) +
  
  # annotations
  annotate("text", x = 2.9, y = -.5,
           label = 'Initial Gradient',
           parse = F, hjust = 0,
           color = 'gray50') +

  annotate("text", x = 3, y = 7,
           label = 'Final Gradient',
           parse = F, hjust = 0,
           angle = 25,
           color = '#000000') +
  
  # add labels to axes
  labs(x = 'Theta',                            
       y = 'y') 

shows the algorithm converging from an initial to final gradient. This simplified version of a gradient descent is fragile and requires manually choosing an \(\alpha\) and threshold for convergence, but demonstrates the idea. The R package numDeriv provides functions to calculate gradients grad(). Of note, our example \(\theta\) was of a single dimension. When the problem includes multiple dimensions of partial second derivatives, they may be collected together in an \(n \times n\) matrix, called a Hessian — the analogue of the second derivative — in R, hessian().

Simplex method

In the above example, we optimized over a continuous variable. Let’s look at another example where the optimal answer must be constrained to integers. For these constraints, we typically use linear programming.

In linear programming, optimization follows a common solution structure. We specify what to minimize or maximize, identify the objective function and constraints (if any, note that the simple least squares problem is a special case of optimization with no constraints). In matrix form,

\[ \text{minimize x}\underbrace{\left[\begin{matrix} c_{1} \\ \vdots \\ c_{n} \end{matrix} \right] ^\mathrm{T} \left[\begin{matrix} x_{1} \\ \vdots \\ x_{n} \end{matrix} \right] }_{\text{Objective}} \text{s.t.}\underbrace{\left[\begin{matrix} a_{11} & \dots & a_{1n} \\ \vdots & \ddots & \vdots \\ a_{m1} & \dots & a_{mn} \end{matrix} \right] \left[\begin{matrix} x_{1} \\ \vdots \\ x_{n} \end{matrix} \right] \geq\left[\begin{matrix} b_{1} \\ \vdots \\ b_{n} \end{matrix} \right]}_{\text{First constraint}},\underbrace{\left[\begin{matrix} x_{1} \\ \vdots \\ x_{n} \end{matrix} \right] \geq 0}_{\text{Second constraint}} \]

Or as restated,

\[ \begin{aligned} \text{min } z &= \boldsymbol{c}\boldsymbol{x} \\ \text{subject to } \boldsymbol{A}\boldsymbol{x} &= \boldsymbol{b} \\ \boldsymbol{x} &\geq 0 \end{aligned} \] The number of functional constraints may be less than the number of variables (i.e., \(n < m\)); thus, we separate \(\boldsymbol{A}=[\boldsymbol{B}, \boldsymbol{N}]\) where \(\boldsymbol{B}\) is a feasible basis. We re-write our problem as,

\[ \begin{aligned} \text{min } \boldsymbol{z} &= \boldsymbol{c}_B \boldsymbol{x}_B + \boldsymbol{c}_N \boldsymbol{x}_N \\ \text{subject to } \boldsymbol{B} \boldsymbol{x}_B + \boldsymbol{N} \boldsymbol{x}_N &= \boldsymbol{b} \\ \boldsymbol{x}_B, \boldsymbol{x}_N &\geq 0 \end{aligned} \]

Solving for \(\boldsymbol{x}_B\), \(\mathbb{I} \boldsymbol{x}_B = \boldsymbol{B}^{-1} \boldsymbol{b} - \boldsymbol{B}^{-1} \boldsymbol{N} \boldsymbol{x}_N\). If \(\boldsymbol{x}_B \geq 0\) (i.e., a feasible solution), we set \(\boldsymbol{x}_N = 0\) (because any other number would increase \(\boldsymbol{z}\)), and substitute both \(\boldsymbol{x}_B\) and \(\boldsymbol{x}_N\) into our objective function, solving for \(\boldsymbol{z}\).

We iterate through each possible combination of bases such that the basis that minimizes \(\boldsymbol{z}\) given the constraints is our optimum. Consider the following example,

\[ \begin{aligned} \text{minimize } z &= 3x_1+x_2+9x_3+x_4 \\ \\ \text{subject to } x_1 + 2x_3 + x_4 &= 4\\ x_2 + x_3 - x_4 &= 2\\ x_i &\geq 0 \end{aligned} \]

In R,

# create a function for our "brute force" simplex method
simplex <- function(A, b, c_x) {
  # number of basis (m) and non-basis (n) variables
  m <- dim(A)[1]
  n <- dim(A)[2] - m
  
  # set up combinations of basis variables
  combs <- t( combn(1:(m + n), n) )
  
  # loop through possible combinations to find lowest cost
  z <- array()
  solution <- array()
  
  for( i in 1:choose(m + n, m) ) {
    
    # get basis variables and related objective constraints
    B <- A[ 1:m, combs[i,] ]
    
    # solve x_B
    x_B <- solve(B) %*% b
    
    # If all x_B >= 0 (feasible solution), find cost
    if( all(x_B >= 0) ) {
      c_B <- c_x[ combs[i,] ]
      
      # get non-basis variables and related objective constraints
      N <- A[ 1:m, setdiff(1:(m + n), combs[i,]) ]
      c_N <- c_x[ setdiff(1:(m + n), combs[i,]) ]
      
      # set non-basis variables to zero
      x_N <- rep(0, n)
      
      # solve for z in objective function
      z[i] <- c_B %*% solve(B) %*% b + 
              ( c_N - c_B %*% solve(B) %*% N ) %*% x_N
      
      # keep minimum cost solution
      if ( z[i] == min(z) ) solution <- x_B 
    }
    else z[i] <- Inf
    
  }
  
  # return NULL if no solution found
  if( min(z) == Inf ) { print("No soluion"); return(NULL) }
  
  # handle multiple solutions
  if( sum( z == min(z) ) > 1 ) { 
    print("Multiple solutions"); 
    return(NULL) 
    }
    
  # return unique solution
  x <- rep(0, m + n)
  x[ combs[which.min(z),] ] <- solution
  s <- list(x = x, 
            z = min(z))
  return(s)
}

# constraint matrix: Ax = b
A <- rbind(c(1, 0, 2, 1),
           c(0, 1, 1, -1))

# right hand side
b <- c(4, 2)

# objective function
c_x <- c(3, 1, 9, 1)

# find an optimum
simplex(A, b, c_x)

## $x
## [1] 0 6 0 4
## 
## $z
## [1] 10

Our above version of a simplex method is brute force. We perform mutiple matrix operations for each possible combination of bases choose(m + n, m) but this gets large quick. A bases with 15 of 30 variables, for example, require the above function to iterate 155,117,520 times! Various packages are available with much more stable and efficient linear programming functions, including lp() in the package lpSolve,

x <- lpSolve::lp(direction = "min", 
                 objective.in = c_x,
                 const.mat = A,
                 const.dir = c("=", "="),
                 const.rhs = b)
x$solution

## [1] 0 6 0 4

An introductory-level review of optimization, beginning with discrete variables and the simplex method is found in Pedregal, Introduction to Optimization (Springer 2004). Optimization over continuous variables is well-presented with modest mathematics in Nocedal & Wright’s Numerical Optimization (Springer 2006). Rigorous treatments are taught in Pinkham, Analysis, Convexity, and Optimization (Draft September 4, 2014); Boyd & Vandenberghe, Convex Optimization (Cambridge University Press 2014); Bertsimas & Tsitsiklis Introduction to Linear Optimization (Athena Scientific 1997); and Sra, et al. Optimization for Machine Learning (MIT Press 2014).