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
.