SML480: Minimizing Cost Functions

nlminb

nlminb lets us find an argument that minimizes a function. We need to provide it with a starting point from which it’ll look for the argument that minimizes the cost function.

f <- function(x){
  x**2 - 2*x + 1   # (x-1)^2
}

g <- function(x){
  (x - 4)^2
}

Let’s find the x’s that minimize f and g:

nlminb(0, f)$par

## [1] 1

nlminb(0, g)$par

## [1] 4

We can use this to minimize the sum of squared errors when predicting log(gdpPercap) from log(lifeExp):

library(gapminder)
quadr.cost <- function(params){
  sum((log(gapminder$lifeExp) - (params[1] + params[2]*log(gapminder$gdpPercap)))^2)
}

We can now find the params that minimize this:

nlminb(c(0,0), quadr.cost)$par

## [1] 2.864177 0.146549

This is the same result as we get with lm:

lm(log(lifeExp) ~ log(gdpPercap), data = gapminder)

## 
## Call:
## lm(formula = log(lifeExp) ~ log(gdpPercap), data = gapminder)
## 
## Coefficients:
##    (Intercept)  log(gdpPercap)  
##         2.8642          0.1465

Let’s make this nicer by making a function that cooks up a quadr.cost for specific x and y vectors:

make.quadr.cost <- function(x, y){
  quadr.cost <- function(params){
    sum((y - (params[1] + params[2]*x))^2)
  }
}

my.lm <- function(x, y){
  cost <- make.quadr.cost(x, y)
  nlminb(c(0, 0), cost)$par
}

We can now run my.lm with any data. Let’s try our example from before again:

my.lm(log(gapminder$gdpPercap), log(gapminder$lifeExp))

## [1] 2.864177 0.146549

Let’s take a closer look at make.quadr.cost: it returns a function, in which the x and y are fixed: all that varies is params.

cost <- make.quadr.cost(log(gapminder$gdpPercap), log(gapminder$lifeExp))
cost

## function(params){
##     sum((y - (params[1] + params[2]*x))^2)
##   }
## <bytecode: 0x7104790>
## <environment: 0x67d5b78>

We can use this function with any params:

cost(c(0, 0.5))

## [1] 363.9332

cost(c(1.5, 0.5))

## [1] 4297.905

A nice thing here is that you can substitute in the arguments to make.quadr.cost to directly make sense of the return value cost.