Posts

Alguien quería un glm forzando determinados coeficientes >0. Una solución 100% bayesiana no era una opción. Hay varias opciones por ahí. Pero me ha sorprendido que la opción esté disponible en glmnet::glmnet: Filosóficamente, es un tanto sorprendente: de alguna manera, glmnet es glm con prioris alrededor del cero. Los límites superiores e inferiores permiten introducir información a priori adicional no necesariamente compatible con la anterior. Desde el punto de vista de la implementación, tiene sentido que estas opciones estén disponibles. glmnet usa coordinate descent como algoritmo de minimización e introducir restricciones en ese tipo de algoritmos es una trivialidad.

La gráfica representa el volumen de la esfera unidad (eje vertical) en el espacio de dimensión x (eje horizontal). Más aquí (de donde procede la gráfica anterior). Moraleja: en dimensiones altas, hay pocos puntos alrededor de uno concreto; o, dicho de otra manera, los puntos están muy alejados entre sí. Por lo que k-vecinos y otros…

Esta entrada debería ser un comentario más en esta otra, pero voy a abusar del privilegio de ser dueño de la plataforma para promocionarla. Voy a decir cosas que son aproximadamente ciertas. Los detalles de la verdad de todo están en la ayuda y el código de sort y sus métodos. En R hay dos métodos de ordenación: shell y radix. El primero es genérico y el segundo es mejor cuando en el vector hay muchos elementos repetidos (p.e., ordenar el censo por provincias). ...

library(microbenchmark) library(ggplot2) a_int <- sample(10:99, 1e6, replace = T) a_char <- paste("P", a_int, sep = "") res <- microbenchmark( sort_int = sort(a_int), sort_char_radix = sort(a_char, method = "radix"), sort_char = sort(a_char), factor_trick = as.character(sort(as.factor(a_char))), times = 50 ) autoplot(res)

Con datos bajados de aquí: library(MicroDatosEs) library(dplyr) library(microbenchmark) library(ggplot2) censo <- censo2010("MicrodatosCP_NV_per_nacional_3VAR.txt") censo_char <- as.data.frame(censo[, c("CPRO", "SEXO", "ECIVIL", "FACTOR")]) censo_factor <- censo_char censo_factor$CPRO <- factor(censo_factor$CPRO) foo <- function(x) x %>% group_by(CPRO) %>% summarise(res = sum((SEXO == "Mujer") * (ECIVIL == "Divorciado") * FACTOR) / sum(FACTOR) * 100) res <- microbenchmark( char = foo(censo_char), factor = foo(censo_factor), times = 10 ) autoplot(res) Da: ¿No es sorprendente? De hecho, plyr es más rápido que dplyr en este caso si no se usan factores. Notas: El hilo de por qué es así en lugar de otra manera se pierde en código escrito en C++. Para otra vida (mía o de otro). Debo agradecer a Diego Castro el intercambio de ideas, código y perplejidades que dieron pie a todo lo de arriba.

Esta entrada es un (otro, que sumar a este o este) recordatorio de que las XI Jornadas de Usuarios de R están en marcha. Y que serán en Madrid, del 14 al 16 de noviembre, etc. Información toda ella que los enlaces anteriores extienden debidamente. (Además hay una tarifa reducida cuyo plazo termina, aviso, muy, muy pronto.)

Esta entrada responde y complementa Malditas proporciones pequeñas I y II_ _trayendo a colación un artículo que ya mencioné en su día y que cuelgo de nuevo: On the Near Impossibility of Measuring the Returns to Advertising. ¡Atención al teorema de la imposibilidad de la Super Bowl! Y el resumen breve: cada vez estamos abocados a medir efectos más y más pequeños. La fruta que cuelga a la altura de la mano ya está en la fragoneta del rumano. Solo nos queda la morralla y cada vez va a costar más separar grano y paja.

Tienes dos variables aleatorias positivamente correlacionadas, $X$ y $Y$ y una muestra de $n$ parejas de ellas $(x_i, y_i)$. La esperanza de $X$, $E(X)$, es conocida y la de $Y$ no. Obviamente, la puedes estimar haciendo $$ E(Y) \sim \frac{1}{n} \sum_i y_i.$$ Sin embargo, la varianza del estimador $$ E(Y) \sim E(X) \frac{\sum y_i}{\sum x_i}$$ es menor. Tengo una explicación de la intuición de por qué eso es cierto en lugar de no serlo. Pero como no sé si es suficientemente buena, dejo que alguien proponga la suya en los comentarios. ...

Era como y se ha convertido en ¡Qué horror! Coda: En otra página de la Wikipedia en la que he caído después por azar he leído la siguiente frase (que por algún motivo encuentro relevante insertar aquí): Los ríos arrastran sedimentos que consiguen colmatar y rellenar de lodo los lagos. Además, la proliferación de ciertas plantas, como el lirio acuático, los obstruye por completo.

Si ejecutas import numpy as np from sklearn.linear_model import LinearRegression n = 1000 X = np.random.rand(n, 2) Y = np.dot(X, np.array([1, 2])) + 1 + np.random.randn(n) / 2 reg = LinearRegression().fit(X, Y) reg.intercept_ reg.coef_ se obtiene más o menos lo esperado. Pero si añades una columna linealmente dependiente, X = np.column_stack((X, 1 * X[:,1])) ocurren cosas de la más calamitosa especie: Y = np.dot(X, np.array([1, 2, 1])) + 1 + np.random.randn(n) / 2 reg = LinearRegression().fit(X, Y) reg.coef_ # array([ 9.89633991e-01, -1.63740303e+14, 1.63740303e+14]) Comentarios: Diríase que la implementación del modelo lineal en scikit-learn no es la que se estudia por doquier (la prima, la inversa, etc.); sospecho que convierte el error cuadrático en una función que depende de los coeficientes y se la pasan a un optimizador (más o menos) genérico. Supongo que la implementación actual pasa todos las pruebas unitarias. Y sospecho, además, que las pruebas unitarias no las ha planteado un estadístico. Así que tal vez los de scikit-learn no saben que tienen problemas de colinealidad; y si alguien se lo ha comentado, igual no han comprendido el issue. Para la predicción, el problema arriba apuntado no es tal. Aun con coeficientes desaforados y siempre que no haya problemas de precisión numérica, tanto da hacer las cosas como todo el mundo o implementando ocurrencias como la anterior. Pero para todo lo demás (p.e., impacto de variables, etc.), la implementación es de traca y no vale para maldita de Dios la cosa. Aunque a estas alturas de siglo, ¿quién en su sano juicio usa el modelo lineal básico? Además, en la práctica, no hay problemas de colinealidad (ni aproximada o, mucho menos, exacta). Hummm… ¿O sí? Mi colega Cañadas ha escrito una entrada en su blog sobre la codificación de variables categóricas donde denuncia cómo en Python las funciones más habituales crean por defecto columnas linealmente dependientes por defecto (por no omitir el primer nivel). O sea, en Python, si no andas con cuidado, acabas con la suela llena de kk de perro.