-
Notifications
You must be signed in to change notification settings - Fork 1
/
Function-factories.qmd
768 lines (538 loc) · 30.3 KB
/
Function-factories.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
# Fábricas de funciones {#sec-function-factories}
\index{function factories}
```{r, include = FALSE}
source("common.R")
```
## Introducción
Una **fábrica de funciones** es una función que hace funciones. Aquí hay un ejemplo muy simple: usamos una fábrica de funciones (`power1()`) para hacer dos funciones secundarias (`square()` y `cube()`):
```{r}
power1 <- function(exp) {
function(x) {
x ^ exp
}
}
square <- power1(2)
cube <- power1(3)
```
No se preocupe si esto aún no tiene sentido, ¡debería tenerlo al final del capítulo!
\index{manufactured functions} \index{functions!manufactured} Llamaré `square()` y `cube()` **funciones fabricadas**, pero este es solo un término para facilitar la comunicación con otros humanos: desde la perspectiva de R, no son diferentes a las funciones creadas de otra manera.
```{r}
square(3)
cube(3)
```
Ya ha aprendido acerca de los componentes individuales que hacen posibles las fábricas de funciones:
- En la @sec-first-class-functions, aprendiste sobre las funciones de primera clase de R. En R, vinculas una función a un nombre de la misma manera que vinculas cualquier objeto a un nombre: con `<-`.
- En la @sec-function-environments, aprendiste que una función captura (encierra) el entorno en el que se crea.
- En la @sec-execution-environments, aprendió que una función crea un nuevo entorno de ejecución cada vez que se ejecuta. Este entorno suele ser efímero, pero aquí se convierte en el entorno envolvente de la función fabricada.
En este capítulo, aprenderá cómo la combinación no obvia de estas tres funciones conduce a la fábrica de funciones. También verá ejemplos de su uso en visualización y estadísticas.
De las tres principales herramientas de programación funcional (funcionales, fábricas de funciones y operadores de funciones), las fábricas de funciones son las menos utilizadas. En general, no tienden a reducir la complejidad general del código, sino que dividen la complejidad en fragmentos más fáciles de digerir. Las fábricas de funciones también son un bloque de construcción importante para los muy útiles operadores de funciones, sobre los cuales aprenderá en el @sec-function-operators.
### Estructura {.unnumbered}
- La @sec-factory-fundamentals comienza el capítulo con una explicación de cómo funcionan las fábricas de funciones, reuniendo ideas del alcance y los entornos. También verá cómo se pueden usar fábricas de funciones para implementar una memoria para funciones, lo que permite que los datos persistan entre llamadas.
- La @sec-graph-fact ilustra el uso de fábricas de funciones con ejemplos de ggplot2. Verá dos ejemplos de cómo funciona ggplot2 con fábricas de funciones proporcionadas por el usuario y un ejemplo de cómo ggplot2 usa una fábrica de funciones internamente.
- La @sec-stat-fact utiliza fábricas de funciones para abordar tres desafíos de las estadísticas: comprender la transformación de Box-Cox, resolver problemas de máxima verosimilitud y dibujar remuestreos de arranque.
- La @sec-functional-factories muestra cómo puede combinar fábricas de funciones y funcionales para generar rápidamente una familia de funciones a partir de datos.
### Requisitos previos {.unnumbered}
Asegúrese de estar familiarizado con el contenido de las Secciones @sec-first-class-functions (funciones de primera clase), @sec-function-environments (el entorno funcional) y @sec-execution-environments (entornos de ejecución) mencionados anteriormente.
Las fábricas de funciones solo necesitan base R. Usaremos un poco de [rlang](https://rlang.r-lib.org) para mirar dentro de ellas más fácilmente, y usaremos [ggplot2](https://%20ggplot2.tidyverse.org) y [scales](https://scales.r-lib.org) para explorar el uso de fábricas de funciones en la visualización.
```{r setup}
library(rlang)
library(ggplot2)
library(scales)
```
## Fundamentos de fábrica {#sec-factory-fundamentals}
La idea clave que hace que las fábricas de funciones funcionen se puede expresar de manera muy concisa:
> El entorno envolvente de la función fabricada es un entorno de ejecución de la fábrica de funciones.
Solo se necesitan unas pocas palabras para expresar estas grandes ideas, pero se necesita mucho más trabajo para entender realmente lo que esto significa. Esta sección te ayudará a juntar las piezas con exploración interactiva y algunos diagramas.
### Entornos
Empecemos echando un vistazo a `square()` y `cube()`:
```{r}
square
cube
```
Es obvio de dónde viene `x`, pero ¿cómo encuentra R el valor asociado con `exp`? La simple impresión de las funciones fabricadas no es reveladora porque los cuerpos son idénticos; los contenidos del entorno envolvente son los factores importantes. Podemos obtener un poco más de información usando `rlang::env_print()`. Eso nos muestra que tenemos dos entornos diferentes (cada uno de los cuales era originalmente un entorno de ejecución de `power1()`). Los entornos tienen el mismo padre, que es el entorno envolvente de `power1()`, el entorno global.
```{r}
env_print(square)
env_print(cube)
```
`env_print()` nos muestra que ambos entornos tienen un enlace a `exp`, pero queremos ver su valor [^function-factories-1]. Podemos hacerlo obteniendo primero el entorno de la función y luego extrayendo los valores:
[^function-factories-1]: Es probable que una versión futura de `env_print()` resuma mejor el contenido, por lo que no necesita este paso.
```{r}
fn_env(square)$exp
fn_env(cube)$exp
```
Esto es lo que hace que las funciones fabricadas se comporten de manera diferente entre sí: los nombres en el entorno adjunto están vinculados a valores diferentes.
### Convenciones de diagrama
También podemos mostrar estas relaciones en un diagrama:
```{r, echo = FALSE, out.width = NULL}
knitr::include_graphics("diagrams/function-factories/power-full.png")
```
Están sucediendo muchas cosas en este diagrama y algunos de los detalles no son tan importantes. Podemos simplificar considerablemente usando dos convenciones:
- Cualquier símbolo flotante libre vive en el entorno global.
- Cualquier entorno sin un padre explícito hereda del entorno global.
```{r, echo = FALSE, out.width = NULL}
knitr::include_graphics("diagrams/function-factories/power-simple.png")
```
Esta vista, que se centra en los entornos, no muestra ningún vínculo directo entre `cube()` y `square()`. Esto se debe a que el enlace se realiza a través del cuerpo de la función, que es idéntico para ambos, pero no se muestra en este diagrama.
Para terminar, veamos el entorno de ejecución de `square(10)`. Cuando `square()` ejecuta `x ^ exp`, encuentra `x` en el entorno de ejecución y `exp` en su entorno adjunto.
```{r}
square(10)
```
```{r, echo = FALSE, out.width = NULL}
knitr::include_graphics("diagrams/function-factories/power-exec.png")
```
### Evaluación forzada {#sec-forcing-evaluation}
```{=tex}
\index{lazy evaluation}
\index{force()}
```
Hay un error sutil en `power1()` causado por una evaluación perezosa. Para ver el problema necesitamos introducir alguna indirección:
```{r}
x <- 2
square <- power1(x)
x <- 3
```
¿Qué debería devolver `square(2)`? Esperarías que devuelva 4:
```{r}
square(2)
```
Desafortunadamente, no es así porque `x` solo se evalúa con pereza cuando se ejecuta `square()`, no cuando se ejecuta `power1()`. En general, este problema surgirá cada vez que cambie un enlace entre llamar a la función de fábrica y llamar a la función fabricada. Es probable que esto suceda rara vez, pero cuando sucede, conducirá a un verdadero error de rascado de cabeza.
Podemos solucionar este problema **forzando** la evaluación con `force()`:
```{r}
power2 <- function(exp) {
force(exp)
function(x) {
x ^ exp
}
}
x <- 2
square <- power2(x)
x <- 3
square(2)
```
Cada vez que cree una fábrica de funciones, asegúrese de que se evalúen todos los argumentos, usando `force()` según sea necesario si el argumento solo lo usa la función fabricada.
### Funciones con estado {#sec-stateful-funs}
\index{<<-}
\index{copy-on-modify!exceptions}
Las fábricas de funciones también le permiten mantener el estado a través de las invocaciones de funciones, lo que generalmente es difícil de hacer debido al principio de nuevo comienzo descrito en la @sec-fresh-start.
Hay dos cosas que lo hacen posible:
- El entorno envolvente de la función fabricada es único y constante.
- R tiene un operador de asignación especial, `<<-`, que modifica los enlaces en el entorno envolvente.
El operador de asignación habitual, `<-`, siempre crea un enlace en el entorno actual. El **operador de superasignación**, `<<-` vuelve a vincular un nombre existente que se encuentra en un entorno principal.
El siguiente ejemplo muestra cómo podemos combinar estas ideas para crear una función que registre cuántas veces ha sido llamada:
```{r}
new_counter <- function() {
i <- 0
function() {
i <<- i + 1
i
}
}
counter_one <- new_counter()
counter_two <- new_counter()
```
```{r, echo = FALSE, out.width = NULL}
knitr::include_graphics("diagrams/function-factories/counter-1.png")
```
Cuando se ejecuta la función fabricada, `i <<- i + 1` modificará `i` en su entorno adjunto. Debido a que las funciones fabricadas tienen entornos envolventes independientes, tienen recuentos independientes:
```{r}
counter_one()
counter_one()
counter_two()
```
```{r, echo = FALSE, out.width = NULL}
knitr::include_graphics("diagrams/function-factories/counter-2.png")
```
Las funciones con estado se utilizan mejor con moderación. Tan pronto como su función comience a administrar el estado de múltiples variables, es mejor cambiar a R6, el tema del @sec-r6.
### Recolección de basura {#sec-factory-pitfalls}
\index{garbage collector!manufactured functions}
Con la mayoría de las funciones, puede confiar en el recolector de basura para limpiar cualquier objeto temporal grande creado dentro de una función. Sin embargo, las funciones fabricadas se aferran al entorno de ejecución, por lo que deberá desvincular explícitamente cualquier objeto temporal grande con `rm()`. Compara los tamaños de `g1()` y `g2()` en el siguiente ejemplo:
```{r}
f1 <- function(n) {
x <- runif(n)
m <- mean(x)
function() m
}
g1 <- f1(1e6)
lobstr::obj_size(g1)
f2 <- function(n) {
x <- runif(n)
m <- mean(x)
rm(x)
function() m
}
g2 <- f2(1e6)
lobstr::obj_size(g2)
```
### Ejercicios
1. La definición de `force()` es simple:
```{r}
force
```
Why is it better to `force(x)` instead of just `x`?
2. Base R contiene dos fábricas de funciones, `approxfun()` y `ecdf()`. Lea su documentación y experimente para descubrir qué hacen las funciones y qué devuelven.
3. Cree una función `pick()` que tome un índice, `i`, como argumento y devuelva una función con un argumento `x` que subjunte `x` con `i`.
```{r, eval = FALSE}
pick(1)(x)
# should be equivalent to
x[[1]]
lapply(mtcars, pick(5))
# should be equivalent to
lapply(mtcars, function(x) x[[5]])
```
4. Cree una función que cree funciones que calculen el i^th^ [momento central](http://en.wikipedia.org/wiki/Central_moment) de un vector numérico. Puedes probarlo ejecutando el siguiente código:
```{r, eval = FALSE}
m1 <- moment(1)
m2 <- moment(2)
x <- runif(100)
stopifnot(all.equal(m1(x), 0))
stopifnot(all.equal(m2(x), var(x) * 99 / 100))
```
5. ¿Qué pasa si no usas un cierre? Haz predicciones, luego verifica con el siguiente código.
```{r}
i <- 0
new_counter2 <- function() {
i <<- i + 1
i
}
```
6. ¿Qué sucede si usa `<-` en lugar de `<<-`? Haz predicciones, luego verifica con el siguiente código.
```{r}
new_counter3 <- function() {
i <- 0
function() {
i <- i + 1
i
}
}
```
## Fábricas gráficas {#sec-graph-fact}
Comenzaremos nuestra exploración de fábricas de funciones útiles con algunos ejemplos de ggplot2.
### Etiquetado
Uno de los objetivos del paquete [scales](http://scales.r-lib.org) es facilitar la personalización de las etiquetas en ggplot2. Proporciona muchas funciones para controlar los detalles finos de ejes y leyendas. Las funciones del formateador[^function-factories-2] son una clase útil de funciones que facilitan el control de la aparición de roturas de ejes. El diseño de estas funciones inicialmente puede parecer un poco extraño: todas devuelven una función, a la que debe llamar para formatear un número.
[^function-factories-2]: Es un desafortunado accidente de la historia que las escalas usen sufijos de función en lugar de prefijos de función. Eso es porque fue escrito antes de que entendiera las ventajas de autocompletar al usar prefijos comunes en lugar de sufijos comunes.
```{r}
y <- c(12345, 123456, 1234567)
comma_format()(y)
number_format(scale = 1e-3, suffix = " K")(y)
```
En otras palabras, la interfaz principal es una fábrica de funciones. A primera vista, esto parece agregar una complejidad adicional por poca ganancia. Pero permite una buena interacción con las escalas de ggplot2, porque aceptan funciones en el argumento `label`:
```{r, out.width = "24%", fig.align="default", fig.width = 2, fig.height = 2, fig.asp = NULL}
df <- data.frame(x = 1, y = y)
core <- ggplot(df, aes(x, y)) +
geom_point() +
scale_x_continuous(breaks = 1, labels = NULL) +
labs(x = NULL, y = NULL)
core
core + scale_y_continuous(
labels = comma_format()
)
core + scale_y_continuous(
labels = number_format(scale = 1e-3, suffix = " K")
)
core + scale_y_continuous(
labels = scientific_format()
)
```
### Contenedores de histograma
Una característica poco conocida de `geom_histogram()` es que el argumento `binwidth` puede ser una función. Esto es particularmente útil porque la función se ejecuta una vez para cada grupo, lo que significa que puede tener diferentes anchos de bin en diferentes facetas, lo que de otro modo no sería posible.
Para ilustrar esta idea y ver dónde podría ser útil el ancho de bin variable, voy a construir un ejemplo en el que un ancho de bin fijo no es muy bueno.
```{r, fig.width = 6, fig.height = 2.5, fig.asp = NULL, out.width = "90%"}
# construir algunos datos de muestra con números muy diferentes en cada celda
sd <- c(1, 5, 15)
n <- 100
df <- data.frame(x = rnorm(3 * n, sd = sd), sd = rep(sd, n))
ggplot(df, aes(x)) +
geom_histogram(binwidth = 2) +
facet_wrap(~ sd, scales = "free_x") +
labs(x = NULL)
```
Aquí cada faceta tiene el mismo número de observaciones, pero la variabilidad es muy diferente. Sería bueno si pudiéramos solicitar que los anchos de intervalo varíen para obtener aproximadamente el mismo número de observaciones en cada intervalo. Una forma de hacerlo es con una fábrica de funciones que ingresa el número deseado de contenedores (`n`) y genera una función que toma un vector numérico y devuelve un ancho de contenedor:
```{r, fig.width = 6, fig.height = 2.5, fig.asp = NULL, out.width = "90%"}
binwidth_bins <- function(n) {
force(n)
function(x) {
(max(x) - min(x)) / n
}
}
ggplot(df, aes(x)) +
geom_histogram(binwidth = binwidth_bins(20)) +
facet_wrap(~ sd, scales = "free_x") +
labs(x = NULL)
```
Podríamos usar este mismo patrón para envolver las funciones base de R que automáticamente encuentran el llamado binwidth óptimo [^function-factories-3], `nclass.Sturges()`, `nclass.scott()` y `nclass .FD()`:
[^function-factories-3]: ggplot2 no expone estas funciones directamente porque no creo que la definición de optimización necesaria para hacer que el problema sea matemáticamente tratable coincida con las necesidades reales de exploración de datos.
```{r, fig.width = 6, fig.height = 2.5, fig.asp = NULL, out.width = "90%"}
base_bins <- function(type) {
fun <- switch(type,
Sturges = nclass.Sturges,
scott = nclass.scott,
FD = nclass.FD,
stop("Unknown type", call. = FALSE)
)
function(x) {
(max(x) - min(x)) / fun(x)
}
}
ggplot(df, aes(x)) +
geom_histogram(binwidth = base_bins("FD")) +
facet_wrap(~ sd, scales = "free_x") +
labs(x = NULL)
```
### `ggsave()`
\index{ggsave()}
Finalmente, quiero mostrar una fábrica de funciones utilizada internamente por ggplot2. `ggplot2:::plot_dev()` es utilizado por `ggsave()` para pasar de una extensión de archivo (por ejemplo, `png`, `jpeg`, etc.) a una función de dispositivo gráfico (por ejemplo, `png()`, `jpeg( )`). El desafío aquí surge porque los dispositivos gráficos básicos tienen algunas inconsistencias menores que debemos corregir:
- La mayoría tiene `filename` como primer argumento, pero algunos tienen `file`.
- El `width` y `height` de los dispositivos gráficos de trama usan unidades de píxeles por defecto, pero los gráficos vectoriales usan pulgadas.
A continuación se muestra una versión levemente simplificada de `plot_dev()`:
```{r}
plot_dev <- function(ext, dpi = 96) {
force(dpi)
switch(ext,
eps = ,
ps = function(path, ...) {
grDevices::postscript(
file = filename, ..., onefile = FALSE,
horizontal = FALSE, paper = "special"
)
},
pdf = function(filename, ...) grDevices::pdf(file = filename, ...),
svg = function(filename, ...) svglite::svglite(file = filename, ...),
emf = ,
wmf = function(...) grDevices::win.metafile(...),
png = function(...) grDevices::png(..., res = dpi, units = "in"),
jpg = ,
jpeg = function(...) grDevices::jpeg(..., res = dpi, units = "in"),
bmp = function(...) grDevices::bmp(..., res = dpi, units = "in"),
tiff = function(...) grDevices::tiff(..., res = dpi, units = "in"),
stop("Unknown graphics extension: ", ext, call. = FALSE)
)
}
plot_dev("pdf")
plot_dev("png")
```
### Ejercicios
1. Comparar y contrastar `ggplot2::label_bquote()` con `scales::number_format()`
## Fábricas estadísticas {#sec-stat-fact}
Los ejemplos más motivadores para las fábricas de funciones provienen de las estadísticas:
- La transformación Box-Cox.
- Remuestreo Bootstrap.
- Estimación de máxima verosimilitud.
Todos estos ejemplos se pueden abordar sin fábricas de funciones, pero creo que las fábricas de funciones son una buena opción para estos problemas y brindan soluciones elegantes. Estos ejemplos requieren algunos antecedentes estadísticos, así que siéntase libre de omitirlos si no tienen mucho sentido para usted.
### La transformación Box-Cox
\index{Box-Cox transformation}
La transformación de Box-Cox (un tipo de [transformación de potencia](https://en.wikipedia.org/wiki/Power_transform)) es una transformación flexible que a menudo se usa para transformar los datos hacia la normalidad. Tiene un solo parámetro, $\lambda$, que controla la fuerza de la transformación. Podríamos expresar la transformación como una función simple de dos argumentos:
```{r}
boxcox1 <- function(x, lambda) {
stopifnot(length(lambda) == 1)
if (lambda == 0) {
log(x)
} else {
(x ^ lambda - 1) / lambda
}
}
```
Pero volver a formular como una fábrica de funciones facilita la exploración de su comportamiento con `stat_function()`:
```{r, out.width = "50%", fig.align="default", fig.width=4}
boxcox2 <- function(lambda) {
if (lambda == 0) {
function(x) log(x)
} else {
function(x) (x ^ lambda - 1) / lambda
}
}
stat_boxcox <- function(lambda) {
stat_function(aes(colour = lambda), fun = boxcox2(lambda), size = 1)
}
ggplot(data.frame(x = c(0, 5)), aes(x)) +
lapply(c(0.5, 1, 1.5), stat_boxcox) +
scale_colour_viridis_c(limits = c(0, 1.5))
# visualmente, log() parece tener sentido como la transformación
# para lambda = 0; a medida que los valores se hacen cada vez más pequeños, la función
# se acerca cada vez más a una transformación de registro
ggplot(data.frame(x = c(0.01, 1)), aes(x)) +
lapply(c(0.5, 0.25, 0.1, 0), stat_boxcox) +
scale_colour_viridis_c(limits = c(0, 1.5))
```
En general, esto le permite usar una transformación de Box-Cox con cualquier función que acepte una función de transformación unaria: no tiene que preocuparse de que esa función proporcione `...` para pasar argumentos adicionales. También creo que la partición de `lambda` y `x` en dos argumentos de función diferentes es natural, ya que `lambda` juega un papel bastante diferente al de `x`.
### Generadores de arranque
\index{bootstrapping}
Las fábricas de funciones son un enfoque útil para el arranque. En lugar de pensar en un solo arranque (¡siempre necesita más de uno!), puede pensar en un **generador** de arranque, una función que produce un nuevo arranque cada vez que se llama:
```{r}
boot_permute <- function(df, var) {
n <- nrow(df)
force(var)
function() {
col <- df[[var]]
col[sample(n, replace = TRUE)]
}
}
boot_mtcars1 <- boot_permute(mtcars, "mpg")
head(boot_mtcars1())
head(boot_mtcars1())
```
La ventaja de una fábrica de funciones es más clara con un bootstrap paramétrico donde primero tenemos que ajustar un modelo. Podemos hacer este paso de configuración una vez, cuando se llama a la fábrica, en lugar de una vez cada vez que generamos el arranque:
```{r}
boot_model <- function(df, formula) {
mod <- lm(formula, data = df)
fitted <- unname(fitted(mod))
resid <- unname(resid(mod))
rm(mod)
function() {
fitted + sample(resid)
}
}
boot_mtcars2 <- boot_model(mtcars, mpg ~ wt)
head(boot_mtcars2())
head(boot_mtcars2())
```
Uso `rm(mod)` porque los objetos del modelo lineal son bastante grandes (incluyen copias completas de la matriz del modelo y los datos de entrada) y quiero mantener la función fabricada lo más pequeña posible.
### Estimación de máxima verosimilitud {#sec-MLE}
\index{maximum likelihood} \index{optimise()} \index{optim()}
El objetivo de la estimación de máxima verosimilitud (MLE) es encontrar los valores de los parámetros para una distribución que hacen que los datos observados sean más probables. Para hacer MLE, comienza con una función de probabilidad. Por ejemplo, tome la distribución de Poisson. Si conocemos $\lambda$, podemos calcular la probabilidad de obtener un vector $\mathbf{x}$ de valores ($x_1$, $x_2$, ..., $x_n$) multiplicando la función de probabilidad de Poisson como sigue:
$P(\lambda, \mathbf{x}) = \prod_{i=1}^{n} \frac{\lambda ^ {x_i} e^{-\lambda}}{x_i!}$
En estadística, casi siempre trabajamos con el registro de esta función. El logaritmo es una transformación monótona que conserva propiedades importantes (es decir, los extremos se encuentran en el mismo lugar), pero tiene ventajas específicas:
- El registro convierte un producto en una suma, con lo que es más fácil trabajar.
- Multiplicar números pequeños produce números aún más pequeños, lo que hace que la aproximación de punto flotante utilizada por una computadora sea menos precisa.
Apliquemos una transformación logarítmica a esta función de probabilidad y simplifiquemos tanto como sea posible:
$\log(P(\lambda, \mathbf{x})) = \sum_{i=1}^{n} \log(\frac{\lambda^{x_i} e^{-\lambda}}{x_i!})$
$\log(P(\lambda, \mathbf{x})) = \sum_{i=1}^{n} \left( x_i \log(\lambda) - \lambda - \log(x_i!) \right)$
$\log(P(\lambda, \mathbf{x})) = \sum*{i=1}^{n} x_i* \log(\lambda) - \sum{i=1}^{n} \lambda - \sum_{i=1}^{n} \log(x_i!)$
$\log(P(\lambda, \mathbf{x})) = \log(\lambda) \sum*{i=1}^{n} x_i - n* \lambda - \sum{i=1}^{n} \log(x_i!)$
Ahora podemos convertir esta función en una función R. La función R es bastante elegante porque R está vectorizado y, dado que es un lenguaje de programación estadístico, R viene con funciones integradas como log-factorial (`lfactorial()`).
```{r}
lprob_poisson <- function(lambda, x) {
n <- length(x)
(log(lambda) * sum(x)) - (n * lambda) - sum(lfactorial(x))
}
```
Considere este vector de observaciones:
```{r}
x1 <- c(41, 30, 31, 38, 29, 24, 30, 29, 31, 38)
```
Podemos usar `lprob_poisson()` para calcular la probabilidad (registrada) de `x1` para diferentes valores de `lambda`.
```{r}
lprob_poisson(10, x1)
lprob_poisson(20, x1)
lprob_poisson(30, x1)
```
Hasta ahora hemos estado pensando en `lambda` como fijo y conocido y la función nos dijo la probabilidad de obtener diferentes valores de `x`. Pero en la vida real, observamos 'x' y es 'lambda' lo que se desconoce. La verosimilitud es la función de probabilidad vista a través de este lente: queremos encontrar la `lambda` que hace que la `x` observada sea la más probable. Es decir, dada `x`, ¿qué valor de `lambda` nos da el valor más alto de `lprob_poisson`()?
En estadística, destacamos este cambio de perspectiva al escribir $f_{\mathbf{x}}(\lambda)$ en lugar de $f(\lambda, \mathbf{x})$. En R, podemos usar una fábrica de funciones. Proporcionamos `x` y generamos una función con un solo parámetro, `lambda`:
```{r}
ll_poisson1 <- function(x) {
n <- length(x)
function(lambda) {
log(lambda) * sum(x) - n * lambda - sum(lfactorial(x))
}
}
```
(No necesitamos `force()` porque `length()` fuerza implícitamente la evaluación de `x`.)
Una cosa buena de este enfoque es que podemos hacer algunos cálculos previos: cualquier término que solo involucre `x` se puede calcular una vez en la fábrica. Esto es útil porque necesitaremos llamar a esta función muchas veces para encontrar la mejor `lambda`.
```{r}
ll_poisson2 <- function(x) {
n <- length(x)
sum_x <- sum(x)
c <- sum(lfactorial(x))
function(lambda) {
log(lambda) * sum_x - n * lambda - c
}
}
```
Ahora podemos usar esta función para encontrar el valor de `lambda` que maximiza la probabilidad (log):
```{r}
ll1 <- ll_poisson2(x1)
ll1(10)
ll1(20)
ll1(30)
```
En lugar de prueba y error, podemos automatizar el proceso de encontrar el mejor valor con `optimise()`. Evaluará `ll1()` muchas veces, utilizando trucos matemáticos para reducir el valor más grande lo más rápido posible. Los resultados nos dicen que el valor más alto es `-30.27` que ocurre cuando `lambda = 32.1`:
```{r}
optimise(ll1, c(0, 100), maximum = TRUE)
```
Ahora, podríamos haber resuelto este problema sin usar una fábrica de funciones porque `optimise()` pasa `...` a la función que se está optimizando. Eso significa que podríamos usar la función de probabilidad de registro directamente:
```{r}
optimise(lprob_poisson, c(0, 100), x = x1, maximum = TRUE)
```
La ventaja de usar una fábrica de funciones aquí es bastante pequeña, pero hay dos sutilezas:
- Podemos precalcular algunos valores en fábrica, ahorrando tiempo de cálculo en cada iteración.
- El diseño de dos niveles refleja mejor la estructura matemática del problema subyacente.
Estas ventajas aumentan en problemas MLE más complejos, donde tiene múltiples parámetros y múltiples vectores de datos.
<!-- GVW: retrocediendo, ¿qué patrones en el código existente deberían buscar las personas que sugieran "Oye, tal vez use una fábrica de funciones aquí"? -->
### Ejercicios
1. En `boot_model()`, ¿por qué no necesito forzar la evaluación de `df` o `formula`?
2. ¿Por qué podrías formular la transformación de Box-Cox de esta manera?
```{r}
boxcox3 <- function(x) {
function(lambda) {
if (lambda == 0) {
log(x)
} else {
(x ^ lambda - 1) / lambda
}
}
}
```
3. ¿Por qué no debe preocuparse de que `boot_permute()` almacene una copia de los datos dentro de la función que genera?
4. ¿Cuánto tiempo ahorra `ll_poisson2()` en comparación con `ll_poisson1()`? Use `bench::mark()` para ver cuánto más rápido ocurre la optimización. ¿Cómo cambia la longitud de 'x' los resultados?
## Fábricas de funciones + funcionales {#sec-functional-factories}
Para terminar el capítulo, mostraré cómo puede combinar funcionales y fábricas de funciones para convertir datos en muchas funciones. El siguiente código crea muchas funciones de potencia con nombres especiales al iterar sobre una lista de argumentos:
```{r}
names <- list(
square = 2,
cube = 3,
root = 1/2,
cuberoot = 1/3,
reciprocal = -1
)
funs <- purrr::map(names, power1)
funs$root(64)
funs$root
```
Esta idea se extiende de manera directa si su fábrica de funciones toma dos (reemplace `map()` con `map2()`) o más (reemplace con `pmap()`) argumentos.
```{=tex}
\index{with()}
\index{attach()}
\index{env\_bind\_*}
```
Una desventaja de la construcción actual es que tienes que prefijar cada llamada de función con `funs$`. Hay tres formas de eliminar esta sintaxis adicional:
- Para un efecto muy temporal, puedes usar `with()`:
```{r}
with(funs, root(100))
```
Recomiendo esto porque deja muy claro cuándo se ejecuta el código en un contexto especial y cuál es ese contexto.
- Para un efecto más prolongado, puede `attach()` las funciones a la ruta de búsqueda, luego `detach()` cuando haya terminado:
```{r}
attach(funs)
root(100)
detach(funs)
```
Probablemente le hayan dicho que evite usar `attach()`, y ese es generalmente un buen consejo. Sin embargo, la situación es un poco diferente a lo habitual porque adjuntamos una lista de funciones, no un data frame. Es menos probable que modifique una función que una columna en un data frame, por lo que algunos de los peores problemas con `attach()` no se aplican.
- Finalmente, podrías copiar las funciones al entorno global con `env_bind()` (aprenderás sobre `!!!` en la @sec-tidy-dots). Esto es en su mayoría permanente:
```{r}
rlang::env_bind(globalenv(), !!!funs)
root(100)
```
Más tarde puede desvincular esos mismos nombres, pero no hay garantía de que no se hayan vuelto a vincular mientras tanto, y es posible que esté eliminando un objeto que otra persona creó.
```{r}
rlang::env_unbind(globalenv(), names(funs))
```
Aprenderá un enfoque alternativo para el mismo problema en la @sec-new-function. En lugar de usar una fábrica de funciones, puede construir la función con cuasicomillas. Esto requiere conocimientos adicionales, pero genera funciones con cuerpos legibles y evita la captura accidental de objetos grandes en el alcance adjunto. Usamos esa idea en la @sec-tag-functions cuando trabajamos en herramientas para generar HTML desde R.
### Ejercicios
1. ¿Cuál de los siguientes comandos es equivalente a `with(x, f(z))`?
(a) `x$f(x$z)`.
(b) `f(x$z)`.
(c) `x$f(z)`.
(d) `f(z)`.
(e) It depends.
2. Compare y contraste los efectos de `env_bind()` frente a `attach()` para el siguiente código.
```{r}
funs <- list(
mean = function(x) mean(x, na.rm = TRUE),
sum = function(x) sum(x, na.rm = TRUE)
)
attach(funs)
mean <- function(x) stop("Hi!")
detach(funs)
env_bind(globalenv(), !!!funs)
mean <- function(x) stop("Hi!")
env_unbind(globalenv(), names(funs))
```