La función “optim” en R

Llega el día en que nos ponen a hacer nuestra primera tarea de econometría: minimizar la suma de los errores (al cuadrado) que se obtienen al ajustar una línea (o un plano) a ciertos datos. Es decir, nuestro querido profesor nos pide minimizar lo siguiente: S=\sum_{i=1}^{n}e_{i}^2. El profesor nos da, por suerte, únicamente datos para dos variables (Y y X). Bueno, ¿qué hacemos? Sabemos varias cosas:

  1. Con sólo una variable independiente, X, sabemos que la suma de los errores S debe ser una función de dos parámetros: \beta_{1} y \beta_{2}. ¿Por qué? Porque los residuales e_{i} son una función de estos dos parámetros (tomando los datos de las observaciones en X y en Y como fijos): e_{i} = y_{i} - (\beta_{1} + \beta_{2}*x_{i}).
  2. Por ende, S = \sum_{i=1}^{n}(y_{i} - (\beta_{1} + \beta_{2}*x_{i}))^2.
  3. Lo que queremos minimizar, entonces, es la parte derecha de la ecuación anterior. Gracias a nuestros amigos Isaac Newton y Gottfried Leibniz, tenemos la maravillosa herramienta del cálculo diferencial para conseguir ese objetivo. Así pues, obtenemos las dos derivadas parciales de la ecuación anterior -una con respecto a \beta_{1} y otra con respecto a \beta_{2} y las igualamos a 0 (para obtener los valores mínimos). Resolvemos las ecuaciones y obtenemos las bonitas betas.

Bueno, todo lo anterior se puede hacer en R sin quebrarse la cabeza gracias a la función “optim” (que también sirve para maximizar funciones -útil, pues, para métodos de Máxima Verosimilitud). ¿Cómo se utiliza? La función (al menos en su “versión básica”) requiere tres argumentos: los valores donde iniciará a buscar los mínimos (argumento par), la función a minimizar (argumento fn), y los datos que se utilizarán para pasar a esa función (argumento data).

Veamos pues un ejemplo. Primero generemos una base de datos con N = 100 y variables Y y X obtenidas aleatoriamente, las Y obtenidas de una distribución uniforme que corre de -30 a 30, y las X de una distribución normal con media 0 y desviación estándar igual a 2 (ponemos un set.seed para que los resultados sean siempre los mismos).

set.seed(12345)
datos = data.frame(y=runif(100, min=-30, max=30), x=rnorm(100, mean=0, sd=2))

Antes de probar optim, veamos cuáles serían los resultados de la regresión utilizando la función lm (que viene empaquetada con R):

lm(formula = datos$y ~ datos$x)
Coefficients:
(Intercept) datos$x
0.7433 -0.6023

Bueno, tenemos que la constante es igual a 0.743 y la pendiente es igual a -0.602. Regresemos a optim.

Primero, necesitamos decirle a optim qué función queremos minimizar (el truco está en saber crear esta función). Llamaré a esta función sum.errores.cuadrados, esta función tendrá dos argumentos (que llamaré justo igual a los argumentos que pide optim: par y data):


#Nota: par[1] y par[2] indica que nos referimos a distintos parámetros; i.e. beta1 y beta2.
#Si queremos más parámetros entonces hay que indicar par[1], par[2], par[3]... etc.
#Por supuesto esto requeriría más variables x, z, w, etc. (una para cada parámetro).
sum.errores.cuadrados = function(data, par) {
 with(data, sum((par[1] + par[2] * x - y)^2))
}

Y, ahora sí, por fin, utilizamos optim (el cual “inserta” a mi función previa, sum.errores.cuadrados, los valores 0 y 0 como valores iniciales para comenzar la “búsqueda” de mínimos):


resultado = optim(par=c(0,0), fn=sum.errores.cuadrados, data=datos)

Y el resultado es el siguiente (vean que, al igual que lo que obtuvimos con la función lm, con optim también obtenemos 0.743 y -0.602):
> resultado

$par
[1] 0.7424558 -0.6023033

$value
[1] 32045.41

$counts
function gradient
43 NA

$convergence
[1] 0

$message
NULL

El resultado de optim nos dice que cuando los argumentos de la función S son 0.743 y -0.602, el valor de la función es 32045.41, el cual debe ser un mínimo. Es decir, 32045.41=f(0.743, -0.062)=\min(S).

Voilà.

Ejemplo adicional: Aquí dejo un ejemplo para los que quieren usar optim para una regresión logit. Los datos de la variable dependiente, Y, ahora son dicotómicos (en el ejemplo de acá abajo los obtuve aleatoriamente usando la funcion rbinom, con el número de “trials” igual a 1 y probabilidad igual a 0.30).


#Generamos los datos
set.seed(12345)
datos = data.frame(y=rbinom(100, 1, .30), x=rnorm(100, 0, 2))

#Creamos la tediosa función a optimizar (una logistic y su log-likelihood, esta última es el logaritmo del producto de las probabilidades para cada observación que se obtienen de la logistic)
#Como queremos maximizar (pues estamos usando máxima verosimilitud), doy un signo negativo a la función, i.e. -log(prod(... etc. etc.
#Este signo negativo es necesario pues, por default, optim busca mínimos (no máximos).
loglike.logit = function(data, par){
 logistic = 1/(1+exp(-par[1]-par[2]*data$x))
 with(data,-log(prod(logistic^y*(1-logistic)^(1-y))))
}

mle.logit = optim(par=c(0,0), fn=loglike.logit, data=datos)
mle.logit

#Comprobamos que optim funciona comparándolo con la función empaquetada de R:
glm(datos$y ~ datos$x, family="binomial")

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: