r

Esto no es práctico, pero sí bonito; bonito, además, de esa forma inasequible a la chusma

Va de muestrear los números $latex 1, \dots, n$ que tienen asignadas probabilidades $latex p_1, \dots, p_n$. Una manera muy impráctica (en R, basta usar sample) y nada intuitiva de hacerlo es recurriendo a la distribución de Gumbel: library(evd) pes <- runif(5) pes <- pes / sum(pes) gammas <- log(pes) + 2 x <- rgumbel(length(pes)) muestra <- which.max(gammas + x) O, en masa, aplicando get_samples <- function(n){ replicate(n, { x <- rgumbel(length(pes)) which.

Una herramienta para el análisis no paramétrico de series temporales

r
Sí, es un ejemplar de mi colección de rarezas estadísticas, técnicas que no entran dentro del currículo estándar pero que pudieran resultar útiles en algún momento, para algún caso particular. Hoy, perfiles matriciales para series temporales, una técnica que sirve esencialmente, para identificar formas que se repiten en series temporales, como Entiendo además que, como consecuencia, también para señalar aquellos ciclos en que se produzcan perfiles anómalos, para su evaluación. Pero dejo que consultéis la información en, por ejemplo, aquí y aquí.

Distribuciones (¿de renta? ¿solo de renta?) a partir de histogramas

En el primer número de la novísima revista Spanish Journal of Statistics aparece un artículo con un título tentador: Recovering income distributions from aggregated data via micro-simulations. Es decir, un artículo que nos puede permitir, por ejemplo, muestrear lo que la AEAT llama rendimientos a partir de lo que publica (aquí): Uno de los métodos de los que sostienen el ignominioso a mí me funciona está basado en el modelo

Más sobre variables instrumentales con R

r
[El título de esta entrada tiene un + delante porque ya escribí sobre el asunto tiempo atrás.] Con la excusa de la reciente publicación del paquete ivreg (para el ajuste de modelos con variables instrumentales, por si el contexto no lo hace evidente), he mirado a ver quién estaba construyendo y ajustando modelos generativos menos triviales que los míos (véase el enlace anterior) para que quede más claro de qué va la cosa.

Contrariamente a lo que creía recordar, "Hot deck" != LOCF

Imputación (que es algo en lo que muy a regañadientes estoy trabajando estos días). Si de verdad tienes que imputar datos en una tabla (y solo en ese caso), solo hay un criterio: construye un modelo para predecir los valores faltantes en función del resto y reemplaza el NA por la su predicción. El modelo puede ser tan tonto como lm(my_col ~ 1, na.rm = T) que resulta en la popular estrategia de reemplazar los NAs por la media del resto de las observaciones.

Misma p, distinto n, luego...

Tres situaciones. La primera: n <- 20 y <- 15 test <- prop.test(y, n, p = .5) test$p.value # [1] 0.04417134 test$conf.int # 0.5058845 0.9040674 La segunda: n <- 200 y <- 115 test <- prop.test(y, n, p = 0.5) test$p.value #[1] 0.04030497 test$conf.int # 0.5032062 0.6438648 Y la tercera: n <- 2000 y <- 1046 test <- prop.test(y, n, p = 0.5) test$p.value #[1] 0.0418688 test$conf.int # 0.5008370 0.5450738 En resumen:

Por supuesto que tengo más variables que observaciones... ¿y?

He intentado replicar los resultados de la entrada de ayer con GAM (vía mgcv) así (véase el enlace anterior para la definición de los datos): library(mgcv) modelo_gam <- gam( y ~ x + s(id, bs = "re"), data = datos, method = "REML", family = "poisson") Y nada: Error in gam(y ~ x + s(id, bs = "re"), data = datos, method = "REML", : Model has more coefficients than data

Sobremuestreando x (y no y)

Construyo unos datos (artificiales, para conocer la verdad): n <- 10000 x1 <- rnorm(n) x2 <- rnorm(n) probs <- -2 + x1 + x2 probs <- 1 / (1 + exp(-probs)) y <- sapply(probs, function(p) rbinom(1, 1, p)) dat <- data.frame(y = y, x1 = x1, x2 = x2) Construyo un modelo de clasificación (logístico, que hoy no hace falta inventar, aunque podría ser cualquier otro): summary(glm(y ~ x1 + x2, data = dat, family = binomial)) #Call: #glm(formula = y ~ x1 + x2, family = binomial, data = dat) # #Deviance Residuals: # Min 1Q Median 3Q Max #-2.

Cuidado con la aleatoriedad "pochola"

Abundo sobre mi entrada del otro día. Usando números aleatorios hirsutos, n <- 200 x <- runif(n) plot(cumsum(x - .5), type = "l") produce mientras que library(randtoolbox) s <- sobol(n, 1, scrambling = 3) plot(cumsum(s - .5), type = "l") genera que tiene un cariz totalmente distinto.

Aleatoriedad hirsuta, aleatoriedad pochola

Contemplando y comparando y se me han venido a la mente los adjetivos hirsuto y pocholo para calificar las respectivas formas de aleatoriedad que representan. La primera es el resultado del habitual n <- 200 x <- runif(n) y <- runif(n) plot(x, y, pch = 16) mientras que la segunda exige el más sofisticado library(randtoolbox) s <- sobol(n, 2, scrambling = 3) x <- s[,1] y <- s[,2] plot(x, y, pch = 16) Se ve que Sobol quería rellenar más armoniosamente el espacio.

De histogramas a distribuciones (usando la de Burr)

Tengo una entrada perpetuamente pendiente que se pospone, entre otras cosas, porque aún no he encontrado una manera satisfactoria para muestrear histogramas. Una de las vías sería dar con (y ajustar) una distribución subyacente que generase unos histogramas similares. Hoy voy a contar un ejemplo de cómo puede fallar tal estrategia. Por un lado he bajado datos de la distribución de renta en España del INE: Por otro, me he dejado convencer temporalmente de que la distribución de Burr podría ser conveniente para modelar la distribución de ingresos de los hogares (Wikipedia dixit!

Optimización estocástica

r
Una de los proyectos en los que estoy trabajando últimamente está relacionado con un problema de optimización no lineal: tengo un modelo (o una familia de modelos) no lineales con una serie de parámetros, unos datos y se trata de lo que no mercería más explicación: encontrar los que minimizan cierta función de error. Tengo implementadas dos vías: La nls, que usa un optimizador numérico genérico para encontrar esos mínimos. (Nótese que uso nls y no nls porque esa función me queda muy corta).