Monthly Archives: January 2014

Bondad de Ajuste en R usando MLE y K-S.

¿Cuántas veces nos hemos encontrado delante de unos datos y nos hubiera gustado saber si siguen algún tipo de distribución concreta? Generalmente se recurre a técnicas de bondad de ajuste en paquetes ya existentes en R, pero lleva mucho tiempo comprobar *una por una*, con estimación de los posibles parámetros por medio. Escribiremos una función que haga todo esto directamente:

La técnica que va a seguir nuestra función va a ser la siguiente:

1. Comprobará el soporte de los datos para saber qué distribuciones puede ajustar.
2. Para cada una de las distribuciones, estima sus parámetros usando máxima verosimilitud
3. Con los parámetros estimados, realiza un test de K-S de los datos para cada distribución.

El paso 3 no es del todo correcto, pues el estadístico de K-S deja de ser válido cuando se estiman parámetros. En cualquier caso, es bastante ilustrativo realizarlo así, y los resultados no son del todo malos.

Ejemplo de aplicación:

– Generamos, digamos, 1250 datos de \(X\sim Weib(20,3)\)

x <- rweibull(1250,20,3)

- Aplicamos la función sobre los datos

> round(goodfit(x),3)
Parameter 1 Parameter 2 p-value
expo 0.342 0.000 0.000
gam 275.506 94.114 0.000
cau 2.961 0.100 0.000
normal 2.927 0.174 0.000
chi 3.867 0.000 0.000
stu 0.916 0.000 0.000
f 130.240 2.689 0.000
unif 2.161 3.322 0.000
weib 20.519 3.005 0.959

Vemos que la estimación de la Weibull es bastante buena, con un p-valor alto. La función puede no funcionar en algunos casos; al fin y al cabo se recurren a técnicas numéricas de optimización y dependiendo de los datos podría no devolver nada por temas de cotas. El código de la función puede encontrarse aquí

Estimación de máxima verosimilitud usando R

No siempre se tienen expresiones explícitas para los estimadores de máxima verosimilitud, de hecho, estos casos suelen ser la orden y no la excepción. En general, se aplican técnicas de optimización numérica sobre la función de densidad propuesta. Para el ejemplo, supondremos que tenemos datos de \(X\sim Ga(p=25,a=8/10)\). Los estimadores de máxima verosimilitud simultáneos no tienen una expresión explícita.

Usamos el paquete stats4.


> library(stats4)

Generamos algunos datos de prueba:


> x <- rgamma(1000,25,8/10)

Escribimos una pequeña función que únicamente dependa de los parámetros y que devuelva \(-\sum_{i=1}^{n}log(f(x))\).


> LL <- function(p,a){ + R <- dgamma(x,p,a) + return(-sum(log(R))) + }

Finalmente, usamos el comando:

mle(LL,start=list(p=1,a=1))

Internamente mle llama a la función optim, que por métodos numéricos se encarga de estimar nuestros parámetros dados un par de valores iniciales para los mismos. La salida puede devolver warnings, causados por la presencia de valores negativos en algún punto para los parámetros. Para evitar esto se pueden establecer cotas para los valores, quedando así:


mle(LL,start=list(p=1,a=1),lower=c(0,0),upper=c(Inf,Inf))

Análisis de sentimiento usando Twitter API

La API de Twitter facilita bastante la extracción de tweets desde cualquier lenguaje de programación. El paquete ‘twitteR’ de R facilita aún más las cosas. Lo primero de todo es registrar nuestra aplicación en la página de desarrolladores de Twitter, para que nos den keys para la aplicación cliente. Una vez hecho esto, pasamos al meollo del asunto. Nos proponemos minar tweets de… por ejemplo, menciones a Mariano Rajoy (@marianorajoy). Para ello, cargamos los paquetes apropiados.


library(plyr)
library(twitteR)

Pasamos a acreditar nuestros credenciales

reqURL <- "https://api.twitter.com/oauth/request_token" accessURL <- "https://api.twitter.com/oauth/access_token" authURL <- "https://api.twitter.com/oauth/authorize" consumerKey <- "xxxxxxxxxxxxxxxxxxx" consumerSecret <- "xxxxxxxxxxxxxxxxxxxxxxxx" twitCred <- OAuthFactory$new(consumerKey=consumerKey, consumerSecret=consumerSecret, requestURL=reqURL, accessURL=accessURL, authURL=authURL) library(RCurl) options(RCurlOptions = list(cainfo = system.file("CurlSSL", "cacert.pem", package = "RCurl"))) twitCred$handshake() registerTwitterOAuth(twitCred)

Pasamos a la parte divertida, cargamos todos los tweets que nos deje la API de menciones a @marianorajoy.


tweets <- searchTwitter('@marianorajoy,',n=2000)

Nos es necesario ahora encontrar un diccionario de palabras positivas y negativas del español. Después de media hora buscando en Google, desistí y opté por traducir otras que encontré del inglés. Puede haber muchos errores por matizaciones, malas traducciones, etc, pero al menos es una primera buena aproximación. Se adjunta el diccionario por si alguien estuviera interesado. Ahora aplicamos una función (cuyo código, bastante sencillo, también se adjunta) que se limita a contar cuántas palabras positivas tiene cada twit y cuántas negativas. Finalmente hace la diferencia entre las anteriores para obtener una puntuación. Como cabe esperar, este algoritmo ni tiene en cuenta intensidad de palabras ni es capaz de identificar sarcasmos (eso para otro día).


texto <- laply(tweets, function(t) {t$getText}) positive <- scan(positive-words.txt, what='character',comment.char=';') negative <- scan(negative-words.txt, what='character',comment.char=';') resultado <- score.sentiment(texto, positive, negative) resultado$score library(ggplot2) qplot(x=resultado$score,geom="histogram")

El resultado, en el cual nuestro Mariano no sale muy bien parado, es el siguiente:

Negativo es malo, positivo bueno.

Sentimiento hacia Rajoy

El código de la función así como los diccionarios usados pueden bajarse de aquí.

Probabilidad de rachas largas

Dado un experimento Bernouilli, nos planteamos cuál es la longitud esperada de la racha más larga. Para ello se recurre a límite de máximo de distribuciones geométricas (una Gumbel). Los resultados vienen bien explicados en un post de John.D.Cook
(enlace)

El documento original al que hace referencia es “The Surprising Predictability of Long Runs” by Mark F. Schilling, Mathematics Magazine 85 (2012), number 2, pages 141–149. Éste puede ser fácilmente encontrado en el siguiente
enlace