l FUNCIONES BÁSICAS DE R

0.28. Funciones round, ceiling, floor y trunc

Existen 4 funciones útiles para modificar u obtener información de un número,
estas funciones son round, ceiling, floor y trunc.
round(x, digits): sirve para redondear un número según los dígitos in-
dicados.
ceiling(x): entrega el mínimo entero mayor o igual que x.
floor(x): entrega el máximo entero menor o igual que x.
trunc(x): entrega la parte entera de un número x.

Ejemplo
Aplique las funciones round, ceiling, floor y trunc a un valor positivo y a
un valor negativo para inspeccionar los resultados.
A continuación el código de prueba para un número positivo cualquiera.

x <- 5.34896 # Número positivo elegido

round(x, digits=3)

## [1] 5.349

ceiling(x)

## [1] 6

floor(x)

## [1] 5

trunc(x)

## [1] 5
A continuación las pruebas con un número negativo cualquiera.

x <- -4.26589 # Número negativo elegido

round(x, digits=3)

## [1] -4.266

ceiling(x)
0.29. FUNCIONES SORT Y RANK li

## [1] -4

floor(x)

## [1] -5

trunc(x)

## [1] -4

0.29. Funciones sort y rank

Las funciones sort y rank son útiles para ordenar los elementos de un vector o
para saber las posiciones que ocuarían los elementos de un vector al ser ordenado.
La estructura de las dos funciones es la siguiente.

sort(x, decreasing = FALSE)

rank(x)

En el parámetro x se ingresa el vector y el parámetro decreasing sirva para

indicar si el ordenamiento es de menor a mayor (por defecto es este) o de mayor
a menor.

Ejemplo
Considere el vector x que tiene los siguientes elementos: 2, 3, 6, 4, 9 y 5. Ordene
el vector de menor a mayor, de mayor a menor y por último encuentre la posición
que ocupan los elementos de x si se ordenaran de menor a mayor.

x <- c(2, 3, 6, 4, 9, 5)
sort(x)

## [1] 2 3 4 5 6 9

sort(x, decreasing=TRUE)

## [1] 9 6 5 4 3 2

rank(x)

## [1] 1 2 5 3 6 4
lii FUNCIONES BÁSICAS DE R

EJERCICIOS
Use funciones o procedimientos (varias líneas) de R para responder cada una de
las siguientes preguntas.

1. ¿Qué cantidad de dinero sobra al repartir 10000$ entre 3 personas?

2. ¿Es el número 4560 divisible por 3?
3. Construya un vector con los números enteros del 2 al 87. ¿Cuáles de esos
números son divisibles por 7?
4. Construya dos vectores, el primero con los números enteros desde 7 hasta
3, el segundo vector con los primeros cinco números positivos divisibles por
5. Sea A la condición de ser par en el primer vector. Sea B la condición
de ser mayor que 10 en el segundo vector. ¿En cuál de las 5 posiciones se
cumple A y B simultáneamente?
5. Consulte este enlace15 en el cual hay una anéctoda de Gauss niño. Use
R para obtener el resultado de la suma solicitada por el profesor del niño
Gauss.
6. Construya un vector con los siguientes elementos: 1, -4, 5, 9, -4. Escriba
un procedimiento para extraer las posiciones donde está el valor mínimo
en el vector.
7. Calcular 8!
𝑖=7
8. Evaluar la siguiente suma ∑𝑖=3 𝑒𝑖
𝑖=10 √
9. Evaluar la siguiente productoria ∏𝑖=1 log 𝑖
10. Construya un vector cualquiera e inviertalo, es decir, que el primer elemen-
to quede de último, el segundo de penúltimo y así sucesivamente. Compare
su resultado con el de la función rev.
11. Create the vector: 1, 2, 3, … , 19, 20.
12. Create the vector: 20, 19, … , 2, 1.
13. Create the vector: 1, −2, 3, −4, 5, −6, … , 19, −20.
14. Create the vector: 0.13 , 0.21 , 0.16 , 0.24 , ..., 0.136 , 0.234 .
100 25 𝑖 𝑖
15. Calculate the following: ∑𝑖=10 (𝑖3 + 4𝑖2 ) and ∑𝑖=1 ( 2𝑖 + 3𝑖2 ).
16. Read the data set available in: [Link]
fhernanb/datos/master/[Link]
17. Use a code to obtain the number of variables of the data set.
18. Use a code to obtain the number of countries in the data set.
19. Which is the country with the higher population?
20. Which is the country with the lowest literacy rate?
21. ¿Qué valor de verdad tiene la siguiente afirmación? “Los resultados de la
función floor y trunc son siempre los mismos”.

En R hay unas bases de datos incluídas, una de ellas es la base de datos llamada
mtcars. Para conocer las variables que están en mtcars usted puede escribir en
la consola ?mtcars o también help(mtcars). De la base mtcars obtenga bases
de datos que cumplan las siguientes condiciones.

15 [Link]
0.29. FUNCIONES SORT Y RANK liii

22. Autos que tengan un rendimiento menor a 18 millas por galón de combus-
tible.
23. Autos que tengan 4 cilindros.
24. Autos que pesen más de 2500 libras y tengan transmisión manual.
liv FUNCIONES BÁSICAS DE R
Instrucciones de control

En R se disponen de varias instrucciones de control para facilitar los procedimien-

tos que un usuario debe realizar. A continuación se explican esas instrucciones
de control.

0.30. Instrucción if
Esta instrucción sirve para realizar un conjunto de operaciones si se cumple
cierta condición. A continuación se muestra la estructura básica de uso.

if (condicion) {
operación 1
operación 2
...
operación final
}

Ejemplo
Una secretaria recibe la información del salario básico semanal de un empleado
y las horas trabajadas durante la semana por ese empleado. El salario básico es
la remuneración por 40 horas de labor por semana, las horas extra son pagadas a
ciencuenta mil pesos. Escriba el procedimiento en R que debe usar la secretaria
para calcular el salario semanal de un empleado que trabajó 45 horas y tiene
salario básico de un millon de pesos.
El código para calcular el salario final del empleado es el siguiente:

sal <- 1 # Salario básico por semana

hlab <- 45 # Horas laboradas por semana

if(hlab > 40) {

hext <- hlab - 40

lv
lvi INSTRUCCIONES DE CONTROL

salext <- hext * 0.05

sal <- sal + salext
}

sal # Salario semanal

## [1] 1.25

0.31. Instrucción if else

Esta instrucción sirve para realizar un conjunto de operaciones cuando NO
se cumple cierta condición evaluada por un if. A continuación se muestra la
estructura básica de uso.

if (condicion) {
operación 1
operación 2
...
operación final
}
else {
operación 1
operación 2
...
operación final
}

0.32. Instrucción ifelse

Se recomienda usar la instrucción ifelse cuando hay una sola instrucción para
el caso if y para el caso else. A continuación se muestra la estructura básica
de uso.

ifelse(condición, operación SI cumple, operación NO cumple)

Ejemplo
Suponga que usted recibe un vector de números enteros, escriba un procedimien-
to que diga si cada elemento del vector es par o impar.
0.33. INSTRUCCIÓN FOR lvii

x <- c(5, 3, 2, 8, -4, 1)

ifelse(x % % 2 == 0, 'Es par', 'Es impar')

## [1] "Es impar" "Es impar" "Es par" "Es par" "Es par" "Es impar"

0.33. Instrucción for

La instrucción for es muy útil para repetir un procedimiento cierta cantidad de
veces. A continuación se muestra la estructura básica de uso.

for (i in secuencia) {
operación 1
operación 2
...
operación final
}

Ejemplo
Escriba un procedimiento para crear 10 muestras de tamaño 100 de una distri-
bución uniforme entre uno y tres. Para cada una de las muestra, se debe contar
el número de elementos de la muestra que fueron mayores o iguales a 2.5.

nrep <- 10 # Número de repeticiones

n <- 100 # Tamaño de la muestra
conteo <- numeric(nrep) # Vector para almacenar el conteo

for (i in 1:nrep) {
x <- runif(n=n, min=1, max=3)
conteo[i] <- sum(x >= 2.5)
}

conteo # Para obtener el conteo

## [1] 24 37 28 26 30 18 29 23 19 19

0.34. Instrucción while

La instrucción while es muy útil para repetir un procedimiento siempre que se
cumple una condición. A continuación se muestra la estructura básica de uso.
lviii INSTRUCCIONES DE CONTROL

while (condición) {
operación 1
operación 2
...
operación final
}

Ejemplo
Suponga que se lanza una moneda en la cual el resultado es cara o sello. Es-
cribir un procedimiento que simule lanzamientos hasta que el número de caras
obtenidas sea 5. El procedimiento debe entregar el historial de lanzamientos.
Para simular el lanzamiento de una moneda se puede usar la función sample
y definiendo el vector resultados con size=1 para simular un lanzamiento, a
continuación el código y tres pruebas ilustrativas.

resultados <- c('Cara', 'Sello')

sample(x=resultados, size=1) # Prueba 1

## [1] "Sello"
Una vez seamos capaces de simular un lanzamiento podemos escribir el procedi-
miento para generar tantos lanzamientos hasta que se cumpla la condición. El
código mostrado abajo permite hacer lo solicitado.

[Link] <- 0 # Contador de lanzamientos

[Link] <- 0 # Contados de caras obtenidas
historial <- NULL # Vector vacío para almacenar

while ([Link] < 5) {

res <- sample(x=resultados, size=1)
[Link] <- [Link] + 1
historial[[Link]] <- res
if (res == 'Cara') {
[Link] <- [Link] + 1
}
}

historial

## [1] "Sello" "Sello" "Sello" "Sello" "Cara" "Cara" "Sello" "Sello" "Cara"
## [10] "Cara" "Cara"
0.35. INSTRUCCIÓN REPEAT lix

[Link]

## [1] 11

La instrucción for se usa cuando sabemos el número de veces que se debe

repetir el procedimiento, mientras que la instrucción while se usa cuando
debemos repetir un procedimiento cuando se cumpla una condición.

0.35. Instrucción repeat

La instrucción while es muy útil para repetir un procedimiento siempre que se
cumple una condición. A continuación se muestra la estructura básica de uso.

repeat {
operación 1
operación 2
...
operación final
if (condición) break
}

Ejemplo
Escribir un procedimiento para ir aumentando de uno en uno el valor de x hasta
que x sea igual a siete El procedimiento debe imprimir por pantalla la secuencia
de valores de x.

x <- 3 # Valor de inicio

repeat {
print(x)
x <- x + 1
if (x == 8) {
break
}
}

## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
lx INSTRUCCIONES DE CONTROL

La instrucción break sirve para salir de un procedimiento iterativo.

Creación de funciones en R

Uno de los atractivos de R es la gran cantidad de funciones que existen para

realizar diversos procedimientos. En este capítulo se explica al lector la forma
de crear sus propias funciones para que pueda realizar diversas tareas y así logre
explotar el potencial que ofrece R.

0.36. ¿Qué es una función en R?

Una función es un conjunto de instrucciones que convierten las entradas (in-
puts) en resultados (outputs) deseados. En la siguiente figura se muestra una
ilustración de lo que es una función.

0.37. Partes de una función en R

Las partes de una función son:
Entradas o argumentos: sirven para ingresar información necesaria para
realizar el procedimiento de la función. Los argumentos pueden estar vacíos
y a la espera de que el usuario ingrese valores, o pueden tener valores por
defecto, esto significa que si el usuario no ingresa un valor, la función usará
el valor por defecto. Una función puede tener o no argumentos de entrada,
en los ejemplos se mostrarán estos casos.
Cuerpo: está formado por un conjunto de instrucciones que transforman
las entradas en las salidas deseadas. Si el cuerpo de la función está formado
por varias instrucciones éstas deben ir entre llaves { }.
Salidas: son los resultados de la función. Toda función debe tener al menos
un resultado. Si una función entrega varios tipos de objetos se acostum-
bra a organizarlos en una lista que puede manejar los diferentes tipos de
objetos.
A continuación se muestra la estructura general de una función en R.

nombre_de_funcion <- function(par1, par2, ...) {

cuerpo

lxi
lxii CREACIÓN DE FUNCIONES EN R

cuerpo
cuerpo
cuerpo
return(resultado)
}

A continuación se mostrarán varios ejemplos sencillos para que el lector aprenda

a construir funciones en R.

Ejemplo
Construir una función que reciba dos números y que entregue la suma de estos
números.
Solución
Lo primero es elegir un nombre apropiado para la función, aquí se usó el nombre
suma porque así se tiene una idea clara de lo que hace la función. La función
suma recibe dos parámetros, x representa el primer valor ingresado mientras que
y representa el segundo. El cuerpo de la función está formado por dos líneas,
en la primera se crea el objeto resultado en el cual se almanacena el valor de
la suma, en la segunda línea se le indica a R que queremos que retorne el valor
de la suma almacenada en el objeto resultado. A continuación se muestra el
código para crear la función solicitada.

suma <- function(x, y) {

resultado <- x + y
return(resultado)
}

Para usar la función creada sólo se debe ejecutar, vamos a obtener la suma de
los valores 4 y 6 usando la función suma, a continuación el código necesario.

suma(x=4, y=6)

## [1] 10
Para funciones simples como la anterior es posible escribirlas en forma más
compacta. Es posible reducir el cuerpo de la función de 2 líneas a sólo una línea
solicitándole a R que retorne directamente la suma sin almacenarla en ningún
objeto. A continuación la función suma modificada.

suma <- function(x, y) {

return(x + y)
}
0.37. PARTES DE UNA FUNCIÓN EN R lxiii

suma(x=4, y=6) # Probando la función

## [1] 10
Debido a que la función suma tiene un cuerpo muy reducido es posible escribirla
en forma más compacta, en una sola línea. A continuación se muestra el código
para reescribir la función.

suma <- function(x, y) x + y

suma(x=4, y=6) # Probando la función

## [1] 10

Ejemplo
Construir una función que genere números aleatorios entre cero y uno hasta
que la suma de éstos números supere por primera vez el valor de 3. La función
debe entregar la cantidad de números aleatorios generados para que se cumpla
la condición.
Solución
Vamos a llamar la función solicitada con el nombre fun1, esta función NO
necesita ningún parámetro de entrada. El valor de 3 que está en la condición
puede ir dentro del cuerpo y por eso no se necesitan parámetros para esta función.
En el cuerpo de la función se genera un vector con un número aleatorio y luego
se chequea si la suma de sus elementos es menor de 3, si se cumple que la suma
es menor que 3 se siguen generando números que se almacenan en el vector num.
Una vez que la suma exceda el valor de 3 NO se ingresa al while y se pide la
longitud del vector o el valor de veces solicitado. A continuación el código de
la función.

fun1 <- function() {

num <- runif(1)
veces <- 1
while (sum(num) < 3) {
veces <- veces + 1
num[veces] <- runif(1)
}
return(veces)
}

fun1() # primera prueba

lxiv CREACIÓN DE FUNCIONES EN R

## [1] 8

Ejemplo
Construir una función que, dado un número entero positivo (cota) ingresado por
el usuario, genere números aleatorios entre cero y uno hasta que la suma de los
números generados exceda por primera vez la cota. La función debe entregar un
vector con los números aleatorios, la suma y la cantidad de números aleatorios.
Si el usuario no ingresa el valor de la cota, se debe asumir igual a 1.

Solución

La función aquí solicitada es similar a la construída en el ejemplo anterior. La

función fun2 tiene un sólo parámetro con el valor por defecto, si el usuario no
ingresa valor a este parámetro, se asumirá el valor de uno. El cuerpo de la función
es similar al anterior. Como la función debe entregar un vector y dos números,
se construye la lista resultado que almacena los tres objetos solicitados. A
continuación el código para función solicitada.

fun2 <- function(cota=1) {

num <- runif(1)
while (sum(num) < cota) {
num <- c(num, runif(1))
}
resultado <- list(vector=num,
suma=sum(num),
cantidad=length(num))
return(resultado)
}

Probando la función con cota de uno.

fun2()

## $vector
## [1] 0.8523376 0.4814579
##
## $suma
## [1] 1.333796
##
## $cantidad
## [1] 2

Probando la función con cota de tres.

0.37. PARTES DE UNA FUNCIÓN EN R lxv

fun2(cota=3)

## $vector
## [1] 0.7703864 0.6567623 0.5173527 0.7785944 0.6926085
##
## $suma
## [1] 3.415704
##
## $cantidad
## [1] 5

Ejemplo
Construya una función que reciba dos números de la recta real y que entregue el
punto médio de estos números. El resultado debe ser un mensaje por pantalla.
Solución
El punto médio entre dos valores es la suma de los números divido entre dos.
La función cat sirve para concatenar objetos y presentarlos por pantalla. A
continuación el código para la función requerida.

medio <- function(a, b) {

medio <- (a + b) / 2
cat("El punto medio de los valores", a, "y", b,
"ingresados es", medio)
}

medio(a=-3, b=-1) # Probando la función

## El punto medio de los valores -3 y -1 ingresados es -2

La función cat es muy útil para presentar resultados por pantalla. Con-
sulte la ayuda de la función para ver otros ejemplos.

EJERCICIOS
Construir funciones en R que realicen lo solicitado.
1. Construya una función que reciba dos números reales a y b, la función
debe decir cuál es el mayor de ellos.
2. Escriba una función llamada media que calcule la media muestral de un
vector numérico x ingresado a la función. A continuación la fórmula para
calcular la media muestral.
lxvi CREACIÓN DE FUNCIONES EN R

𝑛
∑𝑖=1 𝑥𝑖
𝑥̄ =
𝑛

Nota: no puede usar la función mean( ).

3. Construya una función que encuentre las raíces de una ecuación de segundo
grado. El usuario debe suministrar los coeficientes a, b y c de la ecuación
𝑎𝑥2 + 𝑏𝑥 + 𝑐 = 0 y la función debe entregar las raíces.
4. Escribir una función que calcule la velocidad de un proyectil dado que el
usuario ingresa la distancia recorrida en Km y el tiempo necesario en mi-
nutos. Expresar el resultado se debe entregar en metros/segundo, recuerde
que

𝑒𝑠𝑝𝑎𝑐𝑖𝑜
𝑣𝑒𝑙𝑜𝑐𝑖𝑑𝑎𝑑 =
𝑡𝑖𝑒𝑚𝑝𝑜
5. Escribir una función que reciba dos valores 𝑎 y 𝑏 y que los intercambie. Es
decir, si ingresa 𝑎 = 4 y 𝑏 = 9 que la función entregue 𝑎 = 9 y 𝑏 = 4.
6. Construya una función a la cual le ingrese el salario por hora y el número
de horas trabajadas durante una semana por un trabajador. La función
debe calcular el salario neto.
7. Construya una función llamada precio que calcule el precio total de sacar
A fotocopias y B impresiones, sabiendo que los precios son 50 y 100 pesos
para A y B respectivamente si el cliente es un estudiante, y de 75 y 150 para
A y B si el cliente es un profesor. La función debe tener dos argumentos
cuantitativos (A y B) y el argumento lógico estudiante que por defecto
tenga el valor de TRUE. Use la estructura mostrada abajo.

precio <- function(A, B, estudiante=TRUE) {

...
...
...
return([Link])
}

8. Construya una función llamada salario que le ingrese el salario por hora
y el número de horas trabajadas durante una semana por un trabajador.
La función debe calcular el salario neto semanal, teniendo en cuenta que
si el número de horas trabajadas durante la semana es mayor de 48, esas
horas de demás se consideran horas extras y tienen un 35 % de recargo.
Imprima el salario neto. Use la estructura mostrada abajo.

salario <- function([Link], [Link]) {

...
0.37. PARTES DE UNA FUNCIÓN EN R lxvii

...
...
return([Link])
}

9. Construya una función llamada nota que calcule la nota obtenida por un
alumno en una evaluación de tres puntos cuya ponderación o importancia
son 20 %, 30 % y 50 % para los puntos I, II y III respectivamente. Adicio-
nalmente la función debe generar un mensaje sobre si el estudiante aprobó
la evaluación o no. El usuario debe ingresar las notas individuales de los
tres puntos y la función debe entregar la nota final de la evaluación. Use
la estructura mostrada abajo.

nota <- function(p1, p2, p3) {

...
...
...
}

10. Escriba una función llamada minimo que permita obtener el valor mínimo
de un vector numérico. No puede usar ninguna de las funciones básicas de
R como [Link](), [Link](), order(), min( ), max( ), sort( )
u order( ). Use la estructura mostrada abajo.

minimo <- function(x) {

...
...
return(minimo)
}

11. Construya una función que calcule las coordenadas del punto medio 𝑀 en-
tre dos puntos 𝐴 y 𝐵. Vea la siguiente figura para una ilustración. ¿Cuáles
cree usted que deben ser los parámetros de entrada de la función?
lxviii CREACIÓN DE FUNCIONES EN R
Lectura de bases de datos

En este capítulo se mostrará cómo leer una base de datos externa hacia R.

0.38. ¿Qué es una base de datos?

Una base de datos es un arreglo ordenado de variables numéricas, lógicas y
cualitativas.

En la siguiente figura se ilustran los elementos de una base de datos.

0.39. ¿En qué formato almacenar una base de

datos?
Usualmente los archivos con la información para ser leídos por R se pueden
almacenar en formato:

plano con extensión .txt o,

Excel con extensión .csv.

En las secciones siguientes se mostrará cómo almacenar datos en los dos formatos
para ser leídos en R. En el Cuadro 1 se presenta una base de datos pequeña,
tres observaciones y tres variables, que nos servirá como ejemplo para mostrar
cómo se debe almacenar la información.

Cuadro 1: Ejemplo de una base de datos simple.

Edad Fuma Pais

35 TRUE Colombia
46 TRUE Francia
23 FALSE Malta

lxix
lxx LECTURA DE BASES DE DATOS

0.39.1. Almacenamiento de información en Excel

Para almacenar la información del Cuadro 1 en Excel, abrimos un archivo nuevo
archivo de Excel y copiamos la información tal como se muestra en la figura de
abajo. Se debe iniciar en la parte superior izquierda, no se deben dejar filas
vacías, no se debe colorear, no se deben colocar bordes ni nada, se ingresa la
información sin embellecer el contenido. Por último se guarda el archivo en la
carpeta deseada y al momento de nombrar el archivo se debe modificar la opción
tipo de archivo a csv (delimitado por comas).

Recuerde que el archivo de Excel se debe guardar con extensión .csv.

0.39.2. Almacenamiento de información en bloc de notas

Para almacenar la información del Cuadro 1 en un bloc de notas abrimos un
archivo nuevo de bloc de notas y copiamos la información tal como se muestra en
la figura de abajo. Se copian los nombres de las variables o los datos separados
por un espacio obtenido con la tecla tabuladora, cada línea se finaliza con un
enter. Para guardar el archivo se recomienda que el cursor quede al inicio de una
línea vacía. En la figura de abajo se señala la posición del cursor con la flecha
roja, a pesar de que no éxiste línea número 5, el curso debe quedar al inicio de
esa línea número 5.
Es posible mejorar la apariencia de la información almacenada en el bloc de
notas si, en lugar de usar espacios con la barra espaciadora, se colocan los
espacios con la barra tabuladora, así la información se ve más organizada y se
puede chequear fácilmente la información ingresada. En la siguiente figura se
muestra la información para el ejemplo, claramente se nota la organización de
la información.

Una buena práctica es usar la barra tabuladora para separar, eso permite
que la información se vea ordenada.

0.40. Función [Link]

La función [Link] se puede usar para leer bases de datos hacia R. La
estructura de la función con los parámetros más comunes de uso es la siguiente.

[Link](file, header, sep, dec)

Los argumentos de la función [Link] son:

file: nombre o ruta donde están alojados los datos. Puede ser un url o
una dirección del computador. Es también posible usar [Link]()
0.40. FUNCIÓN [Link] lxxi

para que se abra un ventana y adjuntar el archivo deseado manualmente.

header: valor lógico, se usa TRUE si la primera línea de la base de datos
tiene los nombres de las variables, caso contrario se usa FALSE.
sep: tipo de separación interna para los datos dentro del archivo. Los
valores usuales para este parámetros son:
• sep=',' si el archivo tiene extensión .csv.
• sep='' si el archivo es bloc de notas con espacios por la barra espa-
ciadora.
• sep='\t' si el archivo es bloc de notas con espacios por la barra
tabuladora.
dec: símbolo con el cual están indicados los decimales.

Ejemplo
Crear la base de datos del Cuadro 1 en Excel y bloc de notas para practicar la
lectura de base de datos desde R.
Solución
Lo primero que se debe hacer para realizar lo solicitado es construir tres archivos
(uno de Excel y dos bloc de notas) igual a los mostrados en las figuras anterio-
res. Vamos a suponer que los nombres para cada uno de ellos son [Link],
[Link] y [Link] respectivamente.

Para Excel
Para leer el archivo de Excel llamado [Link] podemos usar el siguiente
código.

datos <- [Link](file='C:/Users/mi_usuario/Desktop/[Link]',

header=TRUE, sep=',')
datos

La dirección file='C:/Users/mi_usuario/Desktop/[Link]' le indica a R

en qué lugar del computador debe buscar el archivo, note que se debe usar el
símbolo / para que sea un dirección válida. Substituya la dirección del código
anterior con la dirección donde se encuentra su archivo para que pueda leer la
base de datos.
Si no se conoce la ubicación del archivo a leer o si la dirección es muy extensa,
se puede usar [Link]() para que se abra una ventana y así adjuntar
manualmente el archivo. A continuación se muestra el código para hacerlo de
esta manera.

datos <- [Link]([Link](), header=TRUE, sep=',')

datos
lxxii LECTURA DE BASES DE DATOS

Para bloc de notas con barra espaciadora

Para leer el archivo de Excel llamado [Link] podemos usar el siguiente
código.

datos <- [Link](file='C:/Users/mi_usuario/Desktop/[Link]',

header=TRUE, sep='')
datos

Para bloc de notas con barra tabuladora

Para leer el archivo de Excel llamado [Link] podemos usar el siguiente
código.

datos <- [Link](file='C:/Users/mi_usuario/Desktop/[Link]',

header=TRUE, sep='\t')
datos

El usuario puede usar indiferentemente file='C:/Users/bla/bla' o

[Link]() para ingresar el archivo, con la práctica se aprende a
decidir cuando conviene una u otra forma.

Un error frecuente es escribir la dirección o ubicación del archivo usando

\, lo correcto es usar /.

Ejemplo
Leer la base de datos sobre apartamentos usados en la ciudad de Medellín que
está disponible en la página web cuya url es: [Link]
[Link]/fhernanb/datos/master/aptos2015

Solución

Para leer la base de datos desde una url usamos el siguiente código.

enlace <- '[Link]

datos <- [Link](file=enlace, header=TRUE)

La base de datos ingresada queda en el marco de datos llamado datos y ya está

disponible para usarla.
0.41. LECTURA DE BASES DE DATOS EN EXCEL lxxiii

0.41. Lectura de bases de datos en Excel

Algunas veces los datos están disponibles en un archivo estándar de Excel, y
dentro de cada archivo hojas con la información a utilizar. En estos casos se
recomienda usar el paquete readxl (Wickham and Bryan, 2019) y en particular
la función readxl. A continuación un ejemplo de cómo proceder en estos casos.

Ejemplo
En este enlace16 está disponible un archivo de Excel llamado BD_Excel.xlxs,
una vez se ha abierto la página donde está alojado el archivo, se debe descargar y
guardar en alguna carpeta. El archivo contiene dos bases de datos muy pequeñas,
en la primera hoja llamada Hijos está la información de un grupo de niños y
en la segunda hoja llamada Padres está la información de los padres. ¿Cómo
se pueden leer las dos bases de datos?
Solución
Lo primero que se debe hacer es instalar el paquete readxl, la instalación de
cualquier paquete en un computador se hace una sola vez y éste quedará instala-
do para ser usado las veces que se requiera. La función para instalar un paquete
cualquiera es [Link], a continuación se muestra el código necesario
para instalar el paquete readxl.

[Link]("readxl")

Una vez instalado el paquete es necesario cargarlo, la función para cargar el

paquete en la sesión actual de R es library. La instrucción para cargar el
paquete es la siguiente:

library(readxl)

La instalación de un paquete con [Link] se hace sólo una vez

y no más. Cargar el paquete con library en la sesión actual se debe hacer
siempre que se vaya a usar el paquete.

Luego de haber cargado el paquete readxl se puede usar la función read_xl

para leer la información contenida en las hojas. A continuación el código para
crear la base de datos hijos contenida en el archivo BD_Excel.xlsx.

hijos <- read_excel([Link](), sheet='Hijos')

[Link](hijos) # Para ver el contenido

16 [Link]
lxxiv LECTURA DE BASES DE DATOS

## Edad Grado ComicFav

## 1 8 2 Superman
## 2 6 1 Batman
## 3 9 3 Batman
## 4 10 5 Bob Esponja
## 5 8 4 Batman
## 6 9 4 Bob Esponja
A continuación el código para crear la base de datos padres contenida en el
archivo BD_Excel.xlsx.

padres <- read_excel('BD_Excel.xlsx', sheet='Padres')

[Link](padres) # Para ver el contenido

## Edad EstCivil NumHijos

## 1 45 Soltero 1
## 2 50 Casado 0
## 3 35 Casado 3
## 4 65 Divorciado 1
La función read_excel tiene otros parámetros adicionales útiles para leer bases
de datos, se recomienda consultar la ayuda de la función escribiendo en la consola
help(read_excel).

EJERCICIOS
Realice los siguiente ejercicios propuestos.
1. En el Cuadro 2 se presenta una base de datos sencilla. Almacene la infor-
mación del cuadro en dos archivos diferentes, en Excel y en bloc de notas.
Lea los dos archivos con la función [Link] y compare los resultados
obtenidos con la del Cuadro 2 fuente.
2. En la url [Link]
ter/medidas_cuerpo están disponibles los datos sobre medidas corporales
para un grupo de estudiante de la universidad, use la función [Link]
para leer la base de datos.
0.41. LECTURA DE BASES DE DATOS EN EXCEL lxxv

Cuadro 2: Base de datos para practicar lectura.

Fuma Pasatiempo Num_hermanos Mesada

Si Lectura 0 4500
Si NA 2 2600
No Correr 4 1000
No Correr NA 3990
Si TV 3 2570
No TV 1 2371
Si Correr 1 1389
NA Correr 0 4589
Si Lectura 2 NA
lxxvi LECTURA DE BASES DE DATOS
Tablas de frecuencia

Las tablas de frecuencia son muy utilizadas en estadística y R permite crear ta-
blas de una forma sencilla. En este capítulo se explican las principales funciones
para la elaboración de tablas.

0.42. Tabla de contingencia con table

La función table sirve para construir tablas de frecuencia de una vía, a conti-
nuación la estrctura de la función.

table(..., exclude, useNA)

Los parámetros de la función son:

... espacio para ubicar los nombres de los objetos (variables o vectores)
para los cuales se quiere construir la tabla.
exclude: vector con los niveles a remover de la tabla. Si exclude=NULL
implica que se desean ver los NA, lo que equivale a useNA = 'always'.
useNA: instrucción de lo que se desea con los NA. Hay tres posibles valores
para este parámetro: 'no' si no se desean usar, 'ifany' y 'always' si se
desean incluir.

Ejemplo: tabla de frecuencia de una vía

Considere el vector fuma mostrado a continuación y construya una tabla de
frecuencias absolutas para los niveles de la variable frecuencia de fumar.

fuma <- c('Frecuente', 'Nunca', 'A veces', 'A veces', 'A veces',
'Nunca', 'Frecuente', NA, 'Frecuente', NA, 'hola',
'Nunca', 'Hola', 'Frecuente', 'Nunca')

A continuación se muestra el código para crear la tabla de frecuencias para la

variable fuma.

lxxvii
lxxviii TABLAS DE FRECUENCIA

table(fuma)

## fuma
## A veces Frecuente hola Hola Nunca
## 3 4 1 1 4
De la tabla anterior vemos que NO aparece el conteo de los NA, para obtenerlo
usamos lo siguiente.

table(fuma, useNA='always')

## fuma
## A veces Frecuente hola Hola Nunca <NA>
## 3 4 1 1 4 2
Vemos que hay dos niveles errados en la tabla anterior, Hola y hola. Para
construir la tabla sin esos niveles errados usamos lo siguiente.

table(fuma, exclude=c('Hola', 'hola'))

## fuma
## A veces Frecuente Nunca <NA>
## 3 4 4 2
Por último construyamos la tabla sin los niveles errados y los NA, a esta última
tabla la llamaremos tabla1 para luego poder usarla. Las instrucciones para
hacer esto son las siguientes.

tabla1 <- table(fuma, exclude=c('Hola', 'hola', NA))

tabla1

## fuma
## A veces Frecuente Nunca
## 3 4 4

Al crear una tabla con la instrucción table(var1, var2), la variable 1

quedará por filas mientras que la variable 2 estará en las columnas.

Ejemplo: tabla de frecuencia de dos vías

Considere otro vector sexo mostrado a continuación y construya una tabla de
frecuencias absolutas para ver cómo se relaciona el sexo con fumar del ejemplo
anterior.
0.43. FUNCIÓN [Link] lxxix

sexo <- c('Hombre', 'Hombre', 'Hombre', NA, 'Mujer',

'Casa', 'Mujer', 'Mujer', 'Mujer', 'Hombre', 'Mujer',
'Hombre', NA, 'Mujer', 'Mujer')

Para construir la tabla solicitada usamos el siguiente código.

table(sexo, fuma)

## fuma
## sexo A veces Frecuente hola Hola Nunca
## Casa 0 0 0 0 1
## Hombre 1 1 0 0 2
## Mujer 1 3 1 0 1
De la tabla anterior vemos que aparecen niveles errados en fuma y en sexo, para
retirarlos usamos el siguiente código incluyendo en el parámetro exclude un
vector con los niveles que NO deseamos en la tabla.

tabla2 <- table(sexo, fuma, exclude=c('Hola', 'hola', 'Casa', NA))

tabla2

## fuma
## sexo A veces Frecuente Nunca
## Hombre 1 1 2
## Mujer 1 3 1

0.43. Función [Link]

La función [Link] se utiliza para crear tablas de frecuencia relativa a par-
tir de tablas de frecuencia absoluta, la estructura de la función se muestra a
continuación.

[Link](x, margin=NULL)

x: tabla de frecuencia.
margin: valor de 1 si se desean proporciones por filas, 2 si se desean por
columnas, NULL si se desean frecuencias globales.

Ejemplo: tabla de frecuencia relativa de una vía

Obtener la tabla de frencuencia relativa para la tabla1.
Para obtener la tabla solicitada se usa el siguiente código.
lxxx TABLAS DE FRECUENCIA

[Link](x=tabla1)

## fuma
## A veces Frecuente Nunca
## 0.2727273 0.3636364 0.3636364

Ejemplo: tabla de frecuencia relativa de dos vías

Obtener la tabla de frencuencia relativa para la tabla2.
Si se desea la tabla de frecuencias relativas global se usa el siguiente código. El
resultado se almacena en el objeto tabla3 para ser usado luego.

tabla3 <- [Link](x=tabla2)

tabla3

## fuma
## sexo A veces Frecuente Nunca
## Hombre 0.1111111 0.1111111 0.2222222
## Mujer 0.1111111 0.3333333 0.1111111
Si se desea la tabla de frecuencias relativas marginal por columnas se usa el
siguiente código.

tabla4 <- [Link](x=tabla2, margin=2)

tabla4

## fuma
## sexo A veces Frecuente Nunca
## Hombre 0.5000000 0.2500000 0.6666667
## Mujer 0.5000000 0.7500000 0.3333333

0.44. Función addmargins

Esta función se puede utilizar para agregar los totales por filas o por columnas
a una tabla de frecuencia absoluta o relativa. La estructura de la función es la
siguiente.

addmargins(A, margin)

A: tabla de frecuencia.
margin: valor de 1 si se desean proporciones por columnas, 2 si se desean
por filas, NULL si se desean frecuencias globales.
0.45. FUNCIÓN HIST lxxxi

Ejemplo
Obtener las tablas tabla3 y tabla4 con los totales margines global y por co-
lumnas respectivamente.
Para hacer lo solicitado usamos las siguientes instrucciones.

addmargins(tabla3)

## fuma
## sexo A veces Frecuente Nunca Sum
## Hombre 0.1111111 0.1111111 0.2222222 0.4444444
## Mujer 0.1111111 0.3333333 0.1111111 0.5555556
## Sum 0.2222222 0.4444444 0.3333333 1.0000000

addmargins(tabla4, margin=1)

## fuma
## sexo A veces Frecuente Nunca
## Hombre 0.5000000 0.2500000 0.6666667
## Mujer 0.5000000 0.7500000 0.3333333
## Sum 1.0000000 1.0000000 1.0000000

Note que los valores de 1 y 2 en el parámetro margin de las funciones

[Link] y addmargins significan lo contrario.

0.45. Función hist

Construir tablas de frecuencias para variables cuantitativas es necesario en mu-
chos procedimientos estadísticos, la función hist sirve para obtener este tipo de
tablas. La estructura de la función es la siguiente.

hist(x, breaks='Sturges', [Link]=TRUE, right=TRUE,

plot=FALSE)

Los parámetros de la función son:

x: vector numérico.
breaks: vector con los límites de los intervalos. Si no se especifica se usar
la regla de Sturges para definir el número de intervalos y el ancho.
[Link]: valor lógico, si TRUE una observación 𝑥𝑖 que coincida
con un límite de intervalo será ubicada en el intervalo izquierdo, si FALSE
será incluída en el intervalo a la derecha.
lxxxii TABLAS DE FRECUENCIA

right: valor lógico, si TRUE los intervalos serán cerrados a derecha de la

forma (𝑙𝑖𝑚𝑖𝑛𝑓 , 𝑙𝑖𝑚𝑠𝑢𝑝 ], si es FALSE serán abiertos a derecha.
plot: valor lógico, si FALSE sólo se obtiene la tabla de frecuencias mientras
que con TRUE se obtiene la representación gráfica llamada histograma.

Ejemplo
Genere 200 observaciones aleatorias de una distribución normal con media 𝜇 =
170 y desviación 𝜎 = 5, luego construya una tabla de frecuencias para la muestra
obtenida usando (a) la regla de Sturges y (b) tres intervalos con límites 150, 170,
180 y 190.
Primero se construye el vector x con las observaciones de la distribución normal
por medio de la función rnorm y se especifica la media y desviación solicitada.
Luego se aplica la función hist con el parámetro breaks='Sturges', a conti-
nuación el código utilizado.

x <- rnorm(n=200, mean=170, sd=5)

res1 <- hist(x=x, breaks='Sturges', plot=FALSE)

res1

## $breaks
## [1] 155 160 165 170 175 180 185
##
## $counts
## [1] 4 30 61 71 26 8
##
## $density
## [1] 0.004 0.030 0.061 0.071 0.026 0.008
##
## $mids
## [1] 157.5 162.5 167.5 172.5 177.5 182.5
##
## $xname
## [1] "x"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
El objeto res1 es una lista donde se encuentra la información de la tabla de
frecuencias para x. Esa lista tiene en el elemento breaks los límites inferior y
superior de los intervalos y en el elemento counts están las frecuencias de cada
0.45. FUNCIÓN HIST lxxxiii

uno de los intervalos.

Para obtener las frecuencias de tres intervalos con límites 150, 170, 180 y 190 se
especifica en el parámetros breaks los límites. El código para obtener la segunda
tabla de frecuencias se muestra a continuación.

res2 <- hist(x=x, plot=FALSE,

breaks=c(150, 170, 180, 190))
res2

## $breaks
## [1] 150 170 180 190
##
## $counts
## [1] 95 97 8
##
## $density
## [1] 0.02375 0.04850 0.00400
##
## $mids
## [1] 160 175 185
##
## $xname
## [1] "x"
##
## $equidist
## [1] FALSE
##
## attr(,"class")
## [1] "histogram"

Ejemplo
Construya el vector x con los siguientes elementos: 1.0, 1.2, 1.3, 2.0, 2.5, 2.7,
3.0 y 3.4. Obtenga varias tablas de frecuencia con la función hist variando
los parámetros [Link] y right. Use como límite de los intervalos los
valores 1, 2, 3 y 4.
Lo primero que debemos hacer es crear el vector x solicitado así:

x <- c(1.1, 1.2, 1.3, 2.0, 2.0, 2.5, 2.7, 3.0, 3.4)

En la Figura 1 se muestran los 9 puntos y con color azul se representan los

límites de los intervalos.
A continuación se presenta el código para obtener la tabla de frecuencia usando
lxxxiv TABLAS DE FRECUENCIA

1 2 3 4

Valores de x

Figura 1: Ubicación de los puntos del ejemplo con límites en color azul.

rigth=TRUE, los resultados se almacenan en el objeto res3 y se solicitan sólo

los dos primeros elementos que corresponden a los límites y frecuencias.

res3 <- hist(x, breaks=c(1, 2, 3, 4), right=TRUE, plot=FALSE)

res3[1:2]

## $breaks
## [1] 1 2 3 4
##
## $counts
## [1] 5 3 1
Ahora vamos a repetir la tabla pero usando rigth=FALSE para ver la diferencia,
en res4 están los resultados.

res4 <- hist(x, breaks=c(1, 2, 3, 4), right=FALSE, plot=FALSE)

res4[1:2]

## $breaks
0.45. FUNCIÓN HIST lxxxv

## [1] 1 2 3 4
##
## $counts
## [1] 3 4 2
Al comparar los últimos dos resultados vemos que la primera frecuencia es 5
cuando right=TRUE porque los intervalos se consideran cerrados a la derecha.
Ahora vamos a construir una tabla de frecuencia usando FALSE para los pará-
metros [Link] y right.

res5 <- hist(x, breaks=c(1, 2, 3, 4),

[Link]=FALSE, right=FALSE,
plot=FALSE)
res5[1:2]

## $breaks
## [1] 1 2 3 4
##
## $counts
## [1] 3 4 2
De este último resultado se ve claramente el efecto de los parámetros
[Link] y right en la construcción de tablas de frecuencia.

EJERCICIOS
Use funciones o procedimientos (varias líneas) de R para responder cada una de
las siguientes preguntas.
En el Cuadro 2 se presenta una base de datos sencilla. Lea la base de datos
usando la funcion [Link] y construya lo que se solicita a continuación.
1. Construya una tabla de frecuencia absoluta para la variable pasatiempo.
2. Construya una tabla de frecuencia relativa para la variable fuma.
3. Construya una tabla de frecuencia relativa para las variables pasatiempo
y fuma.
4. ¿Qué porcentaje de de los que no fuman tienen como pasatiempo la lectura.
5. ¿Qué porcentaje de los que corren no fuman?
lxxxvi TABLAS DE FRECUENCIA
Medidas de tendencia
central

En este capítulo se mostrará cómo obtener las diferentes medidas de tendencia

central con R.

Para ilustrar el uso de las funciones se utilizará una base de datos llamada
medidas del cuerpo, esta base de datos cuenta con 6 variables registradas a
un grupo de 36 estudiantes de la universidad. Las variables son:

1. edad del estudiante (años),

2. peso del estudiante (kilogramos),
3. altura del estudiante (centímetros),
4. sexo del estudiante (Hombre, Mujer),
5. muneca: perímetro de la muñeca derecha (centímetros),
6. biceps: perímetro del biceps derecho (centímetros).

A continuación se presenta el código para definir la url donde están los datos,
para cargar la base de datos en R y para mostrar por pantalla un encabezado
(usando head) de la base de datos.

url <- '[Link]

datos <- [Link](file=url, header=T)
head(datos) # Para ver el encabezado de la base de datos

## edad peso altura sexo muneca biceps

## 1 43 87.3 188.0 Hombre 12.2 35.8
## 2 65 80.0 174.0 Hombre 12.0 35.0
## 3 45 82.3 176.5 Hombre 11.2 38.5
## 4 37 73.6 180.3 Hombre 11.2 32.2
## 5 55 74.1 167.6 Hombre 11.8 32.9
## 6 33 85.9 188.0 Hombre 12.4 38.5

lxxxvii
lxxxviii MEDIDAS DE TENDENCIA CENTRAL

0.46. Media
Para calcular la media de una variable cuantitativa se usa la función mean. Los
argumentos básicos de la función mean son dos y se muestran a continuación.

mean(x, [Link] = FALSE)

En el parámetro x se indica la variable de interés para la cual se quiere calcular

la media, el parámetro [Link] es un valor lógico que en caso de ser TRUE, significa
que se deben remover las observaciones con NA, el valor por defecto para este
parámetro es FALSE.

Ejemplo
Suponga que queremos obtener la altura media del grupo de estudiantes.
Para encontrar la media general se usa la función mean sobre el vector númerico
datos$altura.

mean(x=datos$altura)

## [1] 171.5556
Del anterior resultado podemos decir que la estatura media o promedio de los
estudiantes es 171.5555556 centímetros.

Ejemplo
Suponga que ahora queremos la altura media pero diferenciando por sexo.
Para hacer esto se debe primero dividir o partir el vector de altura según los
niveles de la variable sexo, esto se consigue por medio de la función split y el
resultado será una lista con tantos elementos como niveles tenga la variable sexo.
Luego a cada uno de los elementos de la lista se le aplica la función mean con la
ayuda de sapply o tapply. A continuación el código completo para obtener las
alturas medias para hombres y mujeres.

sapply(split(x=datos$altura, f=datos$sexo), mean)

## Hombre Mujer
## 179.0778 164.0333
El resultado es un vector con dos elementos, vemos que la altura media para
hombres es 179.0777778 centímetros y que para las mujeres es de 164.0333333
centímetros.
0.47. MEDIANA lxxxix

¿Qué sucede si se usa tapply en lugar de sapply? Substituya en el código ante-

rior la función sapply por tapply y observe la diferencia entre los resultados.

Ejemplo
Suponga que se tiene el vector edad con las edades de siete personas y supóngase
que para el individuo cinco no se tiene información de su edad, eso significa que
el vector tendrá un NA en la quinta posición.
¿Cuál será la edad promedio del grupo de personas?

edad <- c(18, 23, 26, 32, NA, 32, 29)

mean(x=edad)

## [1] NA
Al correr el código anterior se obtiene un error y es debido al símbolo NA en la
quinta posición. Para calcular la media sólo con los datos de los cuales se tiene
información, se incluye el argumento [Link] = TRUE para que R remueva los NA.
El código correcto a usar en este caso es:

mean(x=edad, [Link]=TRUE)

## [1] 26.66667
De este último resultado se obtiene que la edad promedio de los individuos es
26.67 años.

0.47. Mediana
Para calcular la mediana de una variable cantitativa se usa la función median.
Los argumentos básicos de la función median son dos y se muestran a continua-
ción.

median(x, [Link] = FALSE)

En el parámetro x se indica la variable de interés para la cual se quiere calcular

la mediana, el parámetro [Link] es un valor lógico que en caso de ser TRUE,
significa que se deben remover las observaciones con NA, el valor por defecto
para este parámetro es FALSE.

Ejemplo
Calcular la edad mediana para los estudiantes de la base de datos.
Para obtener la mediana usamos el siguiente código:
xc MEDIDAS DE TENDENCIA CENTRAL

median(x=datos$edad)

## [1] 28
y obtenemos que la mitad de los estudiantes tienen edades mayores o iguales a
28 años.
El resultado anterior se pudo haber obtenido con la función quantile e indi-
cando que se desea el cuantil 50 así:

quantile(x=datos$edad, probs=0.5)

## 50%
## 28

0.48. Moda
La moda de una variable cuantitativa corresponde a valor o valores que más
se repiten, una forma sencilla de encontrar la moda es construir una tabla de
frecuencias y observar los valores con mayor frecuencia.

Ejemplo
Calcular la moda para la variable edad de la base de datos de estudiantes.
Se construye la tabla con la función table y se crea el objeto tabla para alma-
cenarla.

tabla <- table(datos$edad)

tabla

##
## 19 20 21 22 23 24 25 26 28 29 30 32 33 35 37 40 43 45 51 55 65
## 1 1 1 3 2 1 5 3 2 1 2 1 1 2 3 1 2 1 1 1 1
Al mirar con detalle la tabla anterior se observa que el valor que más se repite
es la edad de 25 años en 5 ocasiones. Si la tabla hubiese sido mayor, la ins-
pección visual nos podría tomar unos segundos o hasta minutos y podríamos
equivocarnos, por esa razón es mejor ordenar los resultados de la tabla.
Para observar los valores con mayor frecuencia de la tabla se puede ordenar la
tabla usando la función sort de la siguiente manera:

sort(tabla, decreasing=TRUE)
0.48. MODA xci

##
## 25 22 26 37 23 28 30 35 43 19 20 21 24 29 32 33 40 45 51 55 65
## 5 3 3 3 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1
De esta manera se ve fácilmente que la variable edad es unimodal con valor de
25 años.
xcii MEDIDAS DE TENDENCIA CENTRAL
Medidas de variabilidad

En este capítulo se mostrará cómo obtener las diferentes medidas de variabilidad

con R.
Para ilustrar el uso de las funciones se utilizará la base de datos llamada ap-
tos2015, esta base de datos cuenta con 11 variables registradas a apartamentos
usados en la ciudad de Medellín. Las variables de la base de datos son:
1.
precio: precio de venta del apartamento (millones de pesos),
mt2: área del apartamento (𝑚2 ),
2.
3.
ubicacion: lugar de ubicación del aparamentos en la ciudad (cualitativa),
4.
estrato: nivel socioeconómico donde está el apartamento (2 a 6),
5.
alcobas: número de alcobas del apartamento,
6.
banos: número de baños del apartamento,
7.
balcon: si el apartamento tiene balcón (si o no),
8.
parqueadero: si el apartamento tiene parqueadero (si o no),
9.
administracion: valor mensual del servicio de administración (millones
de pesos),
10. avaluo: valor del apartamento en escrituras (millones de pesos),
11. terminado: si el apartamento se encuentra terminado (si o no).
A continuación se presenta el código para definir la url donde están los datos,
para cargar la base de datos en R y para mostrar por pantalla un encabezado
(usando head) de la base de datos.

url <- '[Link]

datos <- [Link](file=url, header=T)
head(datos) # Para ver el encabezado de la base de datos

## precio mt2 ubicacion estrato alcobas banos balcon parqueadero

## 1 79 43.16 norte 3 3 1 si si
## 2 93 56.92 norte 2 2 1 si si
## 3 100 66.40 norte 3 2 2 no no
## 4 123 61.85 norte 2 3 2 si si
## 5 135 89.80 norte 4 3 2 si no
## 6 140 71.00 norte 3 3 2 no si

xciii
xciv MEDIDAS DE VARIABILIDAD

## administracion avaluo terminado

## 1 0.050 14.92300 no
## 2 0.069 27.00000 si
## 3 0.000 15.73843 no
## 4 0.130 27.00000 no
## 5 0.000 39.56700 si
## 6 0.120 31.14551 si

0.49. Rango
Para calcular el rango de una variable cuantitativa se usa la función range. Los
argumentos básicos de la función range son dos y se muestran abajo.

range(x, [Link] = FALSE)

En el parámetro x se indica la variable de interés para la cual se quiere calcular

el rango, el parámetro [Link] es un valor lógico que en caso de ser TRUE, significa
que se deben remover las observaciones con NA, el valor por defecto para este
parámetro es FALSE.
La función range entrega el valor mínimo y máximo de la variable que se ingresó,
para obtener el valor de rango se debe restar del valor máximo el valor mínimo.

Ejemplo
Suponga que queremos obtener el rango para la variable precio de los aparta-
mentos.
Solución
Para obtener el rango usamos el siguiente código.

range(datos$precio)

## [1] 25 1700
Otra forma de escribir el código anterior de forma secuencial es utilizando el ope-
rador pipe %> %. Este operador se puede leer como “entonces” y permite escribir
código que cuenta una historia.

library(dplyr) # Para cargar el paquete dplyr

datos %> % select(precio) %> % range()

## [1] 25 1700
0.50. VARIANZA xcv

El código de arriba se puede leer así: Tome los datos, entonces seleccione
el precio, entonces calcule el rango.

De la salida anterior podemos ver que los precios de los apartamentos van desde
25 hasta 1700 millones de pesos, es decir, el rango de la variable precio es sería
igual 1700-25=1675 millones de pesos.

Ejemplo
Suponga que queremos obtener nuevamente el rango para la variable precio de
los apartamentos pero diferenciando por el estrato.
Solución
Para calcular el rango (max-min) para el precio pero diferenciando por el estrato
podemos usar el siguiente código.

datos %> %
group_by(estrato) %> %
summarise(el_rango=max(precio)-min(precio))

## # A tibble: 5 x 2
## estrato el_rango
## <int> <dbl>
## 1 2 103
## 2 3 225
## 3 4 610
## 4 5 1325
## 5 6 1560

El código anterior se puede leer así: Tome los datos, entonces agrúpelos
por estrato, entonces haga un resumen llamado el_rango que se obtiene
como el resultado de restar el mínimo de precio al máximo de precio.

De los resultados podemos ver claramente que a medida que aumenta de estra-
to el rango del precio de los apartamentos aumenta. Apartamentos de estrato
bajo tienden a tener precios similares mientras que los precios de venta para
apartamentos de estratos altos tienden a ser muy diferentes entre si.

0.50. Varianza
La varianza es otra medida de qué tanto se alejan las observaciones 𝑥𝑖 en relación
al promedio y se mide en unidades cuadradas. Existen dos formas de calcular
xcvi MEDIDAS DE VARIABILIDAD

la varianza dependiendo de si estamos trabajando con una muestra o con la

población.
La varianza muestral (𝑆 2 ) se define así:
𝑖=𝑛
∑ (𝑥𝑖 − 𝑥)̄ 2
𝑆 = 𝑖=1
2
,
𝑛−1

donde 𝑥̄ representa el promedio muestral.

La varianza poblacional (𝜎2 ) se define así:
𝑖=𝑛
∑𝑖=1 (𝑥𝑖 − 𝜇)2
𝜎2 = ,
𝑛

donde 𝜇 representa el promedio poblacional.

Para calcular la varianza muestral de una variable cuantitativa se usa la función
var. Los argumentos básicos de la función var son dos y se muestran abajo.

var(x, [Link] = FALSE)

En el parámetro x se indica la variable de interés para la cual se quiere calcular

la varianza muestral, el parámetro [Link] es un valor lógico que en caso de
ser TRUE, significa que se deben remover las observaciones con NA, el valor por
defecto para este parámetro es FALSE.

Ejemplo
Suponga que queremos determinar cuál ubicación en la ciudad presenta mayor
varianza en los precios de los apartamentos y cuántos apartamentos hay en cada
ubicación.
Solución
Como nos interesa calcular la varianza y hacer un conteo por cada ubicación,
vamos a agrupar los datos por ubicación. Para realizar lo solicitado podemos
utilizar el siguiente código.

datos %> %
group_by(ubicacion) %> %
summarize(n=n(),
varianza=var(precio))

## # A tibble: 7 x 3
## ubicacion n varianza
## <chr> <int> <dbl>
0.50. VARIANZA xcvii

## 1 aburra sur 169 4169.

## 2 belen guayabal 67 2528.
## 3 centro 38 2588.
## 4 laureles 73 25351.
## 5 norte 10 1009.
## 6 occidente 69 3596.
## 7 poblado 268 84497.
De los resultados anteriores se nota que los apartamentos ubicados en el Poblado
tienen la mayor variabilidad en el precio, este resultado se confirma al dibujar
un boxplot para la variable precio dada la ubicación, en la Figura 2 se muestra el
boxplot y se ve claramente la dispersión de los precios en el Poblado. El código
usado para generar la Figura 2 se presenta a continuación.

with(datos, boxplot(precio ~ ubicacion, ylab='Precio (millones)'))

1500
Precio (millones)

1000
500
0

aburra sur belen guayabal centro laureles norte occidente poblado

ubicacion

Figura 2: Boxplot para el precio de los apartamentos dada la ubicación.

Ejemplo
¿Puedo aplicar la función var y la función sd a marcos de datos?
Solución
La respuesta es NO. La función sd se aplica sólo a vectores mientras que la
función var de puede aplicar tanto a vectores como a marcos de datos. Al ser
aplicada a marcos de datos numéricos se obtiene una matriz en que la diagonal
representa las varianzas de las de cada una de las variables mientras que arriba
y abajo de la diagonal se encuentran las covarianzas entre pares de variables.
Por ejemplo, si aplicamos la función var al marco de datos sólo con las variables
xcviii MEDIDAS DE VARIABILIDAD

precio, área y avaluo se obtiene una matriz de dimensión 3 × 3, a continuación

el código usado.

datos %> %
select(precio, mt2, avaluo) %> %
var()

## precio mt2 avaluo

## precio 61313.15 15874.107 33055.606
## mt2 15874.11 5579.417 9508.188
## avaluo 33055.61 9508.188 28588.853
De la salida anterior se observa que el resultado es una matriz de varianzas y
covarianzas de dimensión 3 × 3.

0.51. Desviación estándar

La desviación estándar es una medida de qué tanto se alejan las observaciones
𝑥𝑖 en relación al promedio y se mide en las mismas unidades de la variable de
interés. Existen dos formas de calcular la desviación estándar dependiendo de
si estamos trabajando con una muestra o con la población.
La desviación estándar muestral (𝑆) para 𝑛 observaciones se define así:

𝑖=𝑛
∑ (𝑥𝑖 − 𝑥)̄ 2
𝑆 = √ 𝑖=1 ,
𝑛−1
donde 𝑥̄ representa el promedio muestral.
La desviación estándar poblacional (𝜎) para 𝑛 observaciones se define así:

𝑖=𝑛
∑𝑖=1 (𝑥𝑖 − 𝜇)2
𝜎=√ ,
𝑛
donde 𝜇 representa el promedio de la población.
Para calcular en R la desviación muestral (𝑆) de una variable cuantitativa se
usa la función sd, los argumentos básicos de la función sd son dos y se muestran
a continuación.

sd(x, [Link] = FALSE)

En el parámetro x se indica la variable de interés para la cual se quiere calcular

la desviación estándar muestral, el parámetro [Link] es un valor lógico que en
caso de ser TRUE, significa que se deben remover las observaciones con NA, el
valor por defecto para este parámetro es FALSE.
0.52. COEFICIENTE DE VARIACIÓN (𝐶𝑉 ) xcix

Ejemplo
Calcular la desviación estándar muestral (𝑆) para la variable precio de los apar-
tamentos.
Solución
Para obtener la desviación estándar muestral (𝑆) solicitada usamos el siguiente
código:

sd(x=datos$precio)

## [1] 247.6149

Ejemplo
Calcular la desviación estándar poblacional (𝜎) para el siguiente conjunto de
5 observaciones: 12, 25, 32, 15, 26.
Solución
Recordemos que las expresiones matemáticas para obtener 𝑆 y 𝜎 son muy si-
milares, la diferencia está en el denominador, para 𝑆 el denominador es 𝑛 − 1
mientras que para 𝜎 es 𝑛. Teniendo esto en cuenta podemos construir nuestra
propia función llamada Sigma que calcule la desviación poblacional. A continua-
ción el código para crear nuestra propia función.

Sigma <- function(x) {

n <- length(x)
desvi <- sqrt(sum((x-mean(x))^2) / n)
return(desvi)
}

Ahora para obtener la desviación estándar poblacional de los datos usamos el

siguiente código.

y <- c(12, 25, 32, 15, 26)

Sigma(y)

## [1] 7.402702

0.52. Coeficiente de variación (𝐶𝑉 )

El coeficiente de variación se define como 𝐶𝑉 = 𝑠/𝑥̄ y es muy sencillo de obte-
nerlo, la función coef_var mostrada abajo permite calcularlo.
c MEDIDAS DE VARIABILIDAD

coef_var <- function(x, [Link] = FALSE) {

sd(x, [Link]=[Link]) / mean(x, [Link]=[Link])
}

Ejemplo
Calcular el 𝐶𝑉 para el vector w definido a continuación.

w <- c(5, -3, NA, 8, 8, 7)

Solución
Vemos que el vector w tiene 6 observaciones y la tercera de ellas es un NA. Lo
correcto aquí es usar la función coef_var definida antes pero indicándole que
remueva los valores faltantes, para eso se usa el siguiente código.

coef_var(x=w, [Link]=T)

## [1] 0.9273618
Medidas de posición

En este capítulo se mostrará cómo obtener las diferentes medidas de posición

con R.
Para ilustrar el uso de las funciones se utilizará una base de datos llamada
medidas del cuerpo, esta base de datos cuenta con 6 variables registradas a
un grupo de 36 estudiantes de la universidad. Las variables son:
1. edad del estudiante (años),
2. peso del estudiante (kilogramos),
3. altura del estudiante (centímetros),
4. sexo del estudiante (Hombre, Mujer),
5. muneca: perímetro de la muñeca derecha (centímetros),
6. biceps: perímetro del biceps derecho (centímetros).
A continuación se presenta el código para definir la url donde están los datos,
para cargar la base de datos en R y para mostrar por pantalla un encabezado
(usando head) de la base de datos.

url <- '[Link]

datos <- [Link](file=url, header=T)
head(datos) # Para ver el encabezado de la base de datos

## edad peso altura sexo muneca biceps

## 1 43 87.3 188.0 Hombre 12.2 35.8
## 2 65 80.0 174.0 Hombre 12.0 35.0
## 3 45 82.3 176.5 Hombre 11.2 38.5
## 4 37 73.6 180.3 Hombre 11.2 32.2
## 5 55 74.1 167.6 Hombre 11.8 32.9
## 6 33 85.9 188.0 Hombre 12.4 38.5

0.53. Cuantiles
Para obtener cualquier cuantil (cuartiles, deciles y percentiles) se usa la función
quantile. Los argumentos básicos de la función quantile son tres y se muestran

ci
cii MEDIDAS DE POSICIÓN

a continuación.

quantile(x, probs, [Link] = FALSE)

En el parámetro x se indica la variable de interés para la cual se quieren calcular

los cuantiles, el parámetro probs sirve para definir los cuantiles de interés y el
parámetro [Link] es un valor lógico que en caso de ser TRUE, significa que se deben
remover las observaciones con NA, el valor por defecto para este parámetro es
FALSE.

Ejemplo
Suponga que queremos obtener el percentil 5, la mediana y el decil 8 pa la altura
del grupo de estudiantes.
Se solicita el percentil 5, la mediana que es el percentil 50 y el decil 8 que
corresponde al percentil 80, por lo tanto es necesario indicarle a la función
quantile que calcule los cuantiles para las ubicaciones 0.05, 0.5 y 0.8, el código
para obtener las tres medidas solicitadas es el siguiente.

quantile(x=datos$altura, probs=c(0.05, 0.5, 0.8))

## 5% 50% 80%
## 155.2 172.7 180.3
Medidas de correlación

En este capítulo se mostrará cómo obtener el coeficiente de correlación lineal

para variables cuantitativas.

0.54. Función cor

La función cor permite calcular el coeficiente de correlación de Pearson, Kendall
o Spearman para dos variables cuantitativas. La estructura de la función es la
siguiente.

cor(x, y, use="everything",
method=c("pearson", "kendall", "spearman"))

Los parámetos de la función son:

x, y: vectores cuantitativos.
use: parámetro que indica lo que se debe hacer cuando se presenten
registros NA en alguno de los vectores. Las diferentes posibilida-
des son: everything, [Link], [Link], [Link] y
[Link], el valor por defecto es everything.
method: tipo de coeficiente de correlación a calcular, por defecto es
pearson, otros valores posibles son kendall y spearman.

Ejemplo
Calcular el coeficiente de correlación de Pearson para las variables área y precio
de la base de datos sobre apartamentos usados.

Lo primero que se debe hacer es cargar la base de datos usando la url apropiada.
Luego de esto se usa la función cor sobre las variables de interés. A continuación
se muestra el código necesario.

ciii
civ MEDIDAS DE CORRELACIÓN

url <- '[Link]

datos <- [Link](file=url, header=T)
cor(x=datos$mt2, y=datos$precio)

## [1] 0.8582585
Del resultado anterior vemos que existe una correlación de 0.8582585 entre las
dos variables, eso significa que apartamentos de mayor área tienden a tener
precios de venta más alto. Este resultado se ilustra en la Figura 3, se nota
claramente que la nube de puntos tiene un pendiente positiva y por eso el signo
del coeficiente de correlación.
A continuación el código para generar la Figura 3.

with(datos, plot(x=mt2, y=precio, pch=20, col='blue',

xlab='Área del apartamento', las=1,
ylab='Precio del apartamento (millones COP)'))

1500
Precio del apartamento (millones COP)

1000

500

100 200 300 400 500

Área del apartamento

Figura 3: Diagrama de dispersión para precio versus área de los apartamentos

usados.

Ejemplo
Para las mismas variables del ejemplo anterior calcular los coeficientes de corre-
lación Kendall y Spearman.
A continuación el código para obtener lo solicitado.
0.54. FUNCIÓN COR cv

cor(x=datos$mt2, y=datos$precio, method='pearson')

## [1] 0.8582585

cor(x=datos$mt2, y=datos$precio, method='kendall')

## [1] 0.6911121

cor(x=datos$mt2, y=datos$precio, method='spearman')

## [1] 0.860306

Ejemplo
Para la base de datos de apartamentos usados, ¿cuáles de las variables cuanti-
tativas tienen mayor correlación?
Lo primero que debemos hacer es determinar cuáles son las cuantitativas de la
base de datos. Para obtener información de las variables que están almacenadas
en el marco de datos llamado datos usamos la función str que muestra la
estructura interna de objeto.

str(datos)

## '[Link]': 694 obs. of 11 variables:

## $ precio : num 79 93 100 123 135 140 145 160 160 175 ...
## $ mt2 : num 43.2 56.9 66.4 61.9 89.8 ...
## $ ubicacion : chr "norte" "norte" "norte" "norte" ...
## $ estrato : int 3 2 3 2 4 3 3 3 4 4 ...
## $ alcobas : int 3 2 2 3 3 3 2 3 4 3 ...
## $ banos : int 1 1 2 2 2 2 2 2 2 2 ...
## $ balcon : chr "si" "si" "no" "si" ...
## $ parqueadero : chr "si" "si" "no" "si" ...
## $ administracion: num 0.05 0.069 0 0.13 0 0.12 0.14 0.127 0 0.123 ...
## $ avaluo : num 14.9 27 15.7 27 39.6 ...
## $ terminado : chr "no" "si" "no" "no" ...
Del anterior resultado vemos que las variables precio, mt2, alcobas, banos, admi-
nistracion y avaluo son las variables cuantitativas, las restantes son cualitativas
(nominal u ordinal). Las posiciones de las variables cuantitativas en el objeto
datos son 1, 2, 5, 6, 9, 10, así podemos construir un marco de datos sólo con la
información cuantitativa, a continuación el código usado.
cvi MEDIDAS DE CORRELACIÓN

[Link] <- datos[, c(1, 2, 5, 6, 9, 10)]

# La siguiente instrucción para editar los nombres de la variables
colnames([Link]) <- c('Precio', 'Área', 'Alcobas',
'Baños', 'Admon', 'Avaluo')
M <- round(cor([Link]), digits=2)
M

## Precio Área Alcobas Baños Admon Avaluo

## Precio 1.00 0.86 0.19 0.63 0.75 0.79
## Área 0.86 1.00 0.31 0.67 0.77 0.75
## Alcobas 0.19 0.31 1.00 0.35 0.16 0.15
## Baños 0.63 0.67 0.35 1.00 0.55 0.53
## Admon 0.75 0.77 0.16 0.55 1.00 0.70
## Avaluo 0.79 0.75 0.15 0.53 0.70 1.00

El anterior resultado representa la matriz de correlaciones entre las variables

cuantitativas, se observa que la mayor correlación es entre las variables precio y
área del apartamento.

Es posible representar gráficamente la matriz de correlaciones M por medio de la

función corrplot del paquete corrplot (Wei and Simko, 2021), a continuación
el código para obtener su representación gráfica.

library('corrplot') # Para cargar el paquete corrplot

## corrplot 0.90 loaded

[Link](M)

En la Figura 4 se muestra la matriz con los coeficientes de correlación. En la

diagonal de la Figura 4 están las variables, por encima están unos círculos de
colores, entre más intensidad del color, ya sea azul o rojo, mayor es la correlación,
colores ténues significan correlación baja; el tamaño de los círculos está asociado
al valor absoluto de correlación. Por debajo de la diagonal se observan los valores
exactos de correlación en colores.

La función corrplot es muy versátil, se pueden obtener diferentes repre-

sentaciones gráficas de la matriz de correlaciones, para conocer las diferen-
tes posibilidades recomendamos consultar este enlace: [Link]
[Link]/web/packages/corrplot/vignettes/corrplot-intro.
html.
0.54. FUNCIÓN COR cvii

Precio 0.8

0.6
0.86 Área
0.4

0.19 0.31 Alcobas 0.2

0.63 0.67 0.35 Baños −0.2

−0.4
0.75 0.77 0.16 0.55 Admon
−0.6

0.79 0.75 0.15 0.53 0.70 Avaluo −0.8

−1

Figura 4: Matriz de coeficientes de correlación.

Ejemplo
Construya dos vectores hipotéticos con el gasto y ahorro de un grupo de 7
familias, incluya dos NA. Calcule el coeficiente de correlación entre ahorro y
gasto, use el parámetro use para manejar los NA.

A continuación se presenta el código para crear los objetos ahorro y gasto con
datos ficticios. Observe que en el primer caso donde se calcula la correlación no
es posible obtener un resultado debido a que por defecto use='everything' y
por lo tanto usa todas las observaciones incluyendo los NA. En el segundo caso si
se obtiene un valor para la correlación debido a que se usó use='[Link]'.

gasto <- c(170, 230, 120, 156, 256, NA, 352)

ahorro <- c(45, 30, NA, 35, 15, 65, 15)

cor(gasto, ahorro)

## [1] NA

cor(gasto, ahorro, use='[Link]')

## [1] -0.8465124
cviii MEDIDAS DE CORRELACIÓN

EJERCICIOS
Use funciones o procedimientos (varias líneas) de R para responder cada una de
las siguientes preguntas.
1. Para cada uno de los estratos socioeconómicos, calcular el coeficiente de
correlación lineal de Pearson para las variables precio y área de la base de
datos de los apartamentos usados.
2. Calcular los coeficientes de correlación Pearson, Kendall y Spearman para
las variables cuantitativas de la base de datos sobre medidas del cuerpo
explicada en el Capítulo 0.45. La url con la información es la siguiente:
[Link]
das_cuerpo
3. Represente gráficamente las matrices de correlación obtenidas en el ejerci-
cio anterior.
Distribuciones discretas

En este capítulo se mostrarán las funciones de R para distribuciones discretas.

0.55. Funciones disponibles para distribuciones

discretas
Para cada distribución discreta se tienen 4 funciones, a continuación el listado
de funciones y su utilidad.

dxxx(x, ...) # Función de masa de probabilidad, f(x)

pxxx(q, ...) # Función de distribución acumulada hasta q, F(x)
qxxx(p, ...) # Cuantil para el cual P(X <= q) = p
rxxx(n, ...) # Generador de números aleatorios.

En el lugar de las letras xxx se de debe colocar el nombre de la distribución

en R, a continuación el listado de nombres disponibles para las 6 distribuciones
discretas básicas.

binom # Binomial
geo # Geométrica
nbinom # Binomial negativa
hyper # Hipergeométrica
pois # Poisson
multinom # Multinomial

Combinando las funciones y los nombres se tiene un total de 24 funciones, por

ejemplo, para obtener la función de masa de probabilidad 𝑓(𝑥) de una binomial
se usa la función dbinom( ) y para obtener la función acumulada 𝐹 (𝑥) de una
Poisson se usa la función ppois( ).

cix
cx DISTRIBUCIONES DISCRETAS

Ejemplo binomial
Suponga que un grupo de agentes de tránsito sale a una vía principal para
revisar el estado de los buses de transporte intermunicipal. De datos históricos
se sabe que un 10 % de los buses generan una mayor cantidad de humo de la
permitida. En cada jornada los agentes revisan siempre 18 buses, asuma que el
estado de un bus es independiente del estado de los otros buses.
1) Calcular la probabilidad de que se encuentren exactamente 2 buses que
generan una mayor cantidad de humo de la permitida.
Aquí se tiene una distribucion 𝐵𝑖𝑛𝑜𝑚𝑖𝑎𝑙(𝑛 = 18, 𝑝 = 0.1) y se desea calcular
𝑃 (𝑋 = 2). Para obtener esta probabilidad se usa la siguiente instrucción.

dbinom(x=2, size=18, prob=0.10)

## [1] 0.2835121
Así 𝑃 (𝑋 = 2) = 0.2835.
2) Calcular la probabilidad de que el número de buses que sobrepasan el
límite de generación de gases sea al menos 4.
En este caso interesa calcular 𝑃 (𝑋 ≥ 4), para obtener esta probabilidad se usa
la siguiente instrucción.

sum(dbinom(x=4:18, size=18, prob=0.10))

## [1] 0.09819684
Así 𝑃 (𝑋 ≥ 4) = 0.0982
3) Calcular la probabilidad de que tres o menos buses emitan gases por enci-
ma de lo permitido en la norma.
En este caso interesa 𝑃 (𝑋 ≤ 3) lo cual es 𝐹 (𝑥 = 3), por lo tanto, la instrucción
para obtener esta probabilidad es

pbinom(q=3, size=18, prob=0.10)

## [1] 0.9018032
Así 𝑃 (𝑋 ≤ 3) = 𝐹 (𝑥 = 3) = 0.9018
4) Dibujar la función de masa de probabilidad.
Para dibujar la función de masa de probabilidad para una 𝐵𝑖𝑛𝑜𝑚𝑖𝑎𝑙(𝑛 = 18, 𝑝 =
0.1) se usa el siguiente código.
0.55. FUNCIONES DISPONIBLES PARA DISTRIBUCIONES DISCRETAS cxi

x <- 0:18 # Soporte (dominio) de la variable

Probabilidad <- dbinom(x=x, size=18, prob=0.1)
plot(x=x, y=Probabilidad,
type='h', las=1, lwd=6)

0.30

0.25

0.20
Probabilidad

0.15

0.10

0.05

0.00

0 5 10 15

Figura 5: Función de masa de probabilidad para una 𝐵𝑖𝑛𝑜𝑚𝑖𝑎𝑙(𝑛 = 18, 𝑝 =

0.1).

En la Figura 5 se muestra la función de masa de probabilidad para la

𝐵𝑖𝑛𝑜𝑚𝑖𝑎𝑙(𝑛 = 18, 𝑝 = 0.1), de esta figura se observa claramente que la mayor
parte de la probabilidad está concentrada para valores pequeños de 𝑋 debido
a que la probabilidad de éxito individual es 𝑝 = 0.10. Valores de 𝑋 ≥ 7 tienen
una probabilidad muy pequeña y es por eso que las longitudes de sus barras
son muy cortas.

5) Generar con 100 de una distribución 𝐵𝑖𝑛𝑜𝑚𝑖𝑎𝑙(𝑛 = 18, 𝑝 = 0.1) y luego

calcular las frecuencias muestrales y compararlas con las probabilidades
teóricas.

La muestra aleatoria se obtiene con la función rbinom y los resultados se alma-

cenan en el objeto m, por último se construye la tabla de frecuencias relativas, a
continuación el código usado.
cxii DISTRIBUCIONES DISCRETAS

m <- rbinom(n=100, size=18, prob=0.1)

m # Para ver lo que hay dentro de m

## [1] 0 1 4 1 3 4 2 3 1 0 3 2 1 4 1 1 1 3 0 5 0 3 5 3 2 3 1 1 2 4 1 3 2 2 3 2 2
## [38] 3 2 3 0 2 2 2 1 1 3 1 1 1 1 1 1 1 1 3 1 2 4 4 7 2 4 0 1 2 1 1 2 3 2 1 3 1
## [75] 1 2 3 5 1 5 0 1 2 1 4 1 1 2 1 2 0 2 2 3 3 0 1 1 3 2

[Link](table(m)) # Tabla de frecuencia relativa

## m
## 0 1 2 3 4 5 7
## 0.09 0.35 0.24 0.19 0.08 0.04 0.01
A pesar de ser una muestra aleatoria de sólo 100 observaciones, se observa que
las frecuencias relativas obtenidas son muy cercanas a las mostradas en la Figura
5.

Ejemplo geométrica
En una línea de producción de bombillos se sabe que sólo el 1 % de los bombillos
son defectuosos. Una máquina automática toma un bombillo y lo prueba, si el
bombillo enciende, se siguen probando los bombillos hasta que se encuentre un
bombillo defectuoso, ahí se para la línea de producción y se toman los correctivos
necesarios para mejorar el proceso.
1) Calcular la probabilidad de que se necesiten probar 125 bombillos para
encontrar el primer bombillo defectuoso.
En la distribución geométrica, la variable 𝑋 representa el número de fracasos
antes de encontrar el único éxito, por lo tanto, en este caso el interés es calcular
𝑃 (𝑋 = 124). La instrucción para obtener esta probabiliad es la siguiente.

dgeom(x=124, prob=0.01)

## [1] 0.002875836
2) Calcular 𝑃 (𝑋 ≤ 8).
En este caso interesa 𝑃 (𝑋 ≤ 50) lo que equivale a 𝐹 (8), la instrucción para
obtener la probabilidad es la siguiente.

pgeom(q=50, prob=0.01)

## [1] 0.401044
3) Encontrar el cuantil 𝑞 tal que 𝑃 (𝑋 ≤ 𝑞) = 0.40.
0.55. FUNCIONES DISPONIBLES PARA DISTRIBUCIONES DISCRETAS cxiii

En este caso interesa encontrar el cuantil 𝑞 que cumpla la condición de que hasta
𝑞 esté el 40 % de las observaciones, por esa razón se usa la función qgeom como
se muestra a continuación.

qgeom(p=0.4, prob=0.01)

## [1] 50

Note que las funciones pxxx y qxxx están relacionadas, pxxx entrega la
probabilidad hasta el cuantil 𝑞 mientras qxxx entrega el cuantil en el que
se acumula 𝑝 probabilidad.

Ejemplo binomial negativa

Una familia desea tener hijos hasta conseguir 2 niñas, la probabilidad individual
de obtener una niña es 0.5 y se supone que todos los nacimientos son individuales,
es decir, un sólo bebé.
1) Calcular la probabilidad de que se necesiten 4 hijos, es decir, 4 nacimientos
para consguir las dos niñas.
En este problema se tiene una distribución binomial negativa con 𝑟 = 2 niñas,
los éxitos deseados por la familia. La variable 𝑋 representa los fracasos, es decir
los niños, hasta que se obtienen los éxitos 𝑟 = 2 deseados.
En este caso lo que interesa es 𝑃 (familia tenga 4), en otras palabras interesa
𝑃 (𝑋 = 2), la instrucción para calcular la probabilidad es la siguiente.

dnbinom(x=2, size=2, prob=0.5)

## [1] 0.1875
2) Calcular 𝑃 (familia tenga al menos 4 hijos).
Aquí interesa calcular 𝑃 (𝑋 ≥ 2) = 𝑃 (𝑋 = 2) + 𝑃 (𝑋 = 3) + …, como esta
probabilidad va hasta infinito, se debe usar el complemento así:

𝑃 (𝑋 ≥ 2) = 1 − [𝑃 (𝑋 = 0) + 𝑃 (𝑋 = 1)]
y para obtener la probabilidad solicitada se puede usar la función dnbinom de
la siguiente manera.

1 - sum(dnbinom(x=0:1, size=2, prob=0.5))

## [1] 0.5
cxiv DISTRIBUCIONES DISCRETAS

Otra forma para obtener la probabilidad solicitada es por medio de la función

pnbinom de la siguiente manera.

1 - pnbinom(q=1, size=2, prob=0.5)

## [1] 0.5

Ejemplo hipergeométrica
Un lote de partes para ensamblar en una empresa está formado por 100 elemen-
tos del proveedor A y 200 elementos del proveedor B. Se selecciona una muestra
de 4 partes al azar sin reemplazo de las 300 para una revisión de calidad.
1) Calcular la probabilidad de que todas las 4 partes de la muestra sean del
proveedor A.
Aquí se tiene una situación que se puede modelar por medio de una distribución
hipergeométrica con 𝑚 = 100 éxitos en la población, 𝑛 = 200 fracasos en la
población y 𝑘 = 4 el tamaño de la muestra. El objetivo es calcular 𝑃 (𝑋 = 4),
para obtener esta probabilidad se usa la siguiente instrucción.

dhyper(x=4, m=100, k=4, n=200)

## [1] 0.01185408
2) Calcular la probabilidad de que dos o más de las partes sean del proveedor
A.
Aquí interesa 𝑃 (𝑋 ≥ 2), la instrucción para obtener esta probabilidad es.

sum(dhyper(x=2:4, m=100, k=4, n=200))

## [1] 0.4074057

Ejemplo Poisson
En una editorial se asume que todo libro de 250 páginas tiene en promedio 50
errores.
1) Encuentre la probabilidad de que en una página cualquiera no se encuen-
tren errores.
Este es un problema de distribución Poisson con tasa promedio de éxitos dada
por:

50 𝑒𝑟𝑟𝑜𝑟𝑒𝑠 0.2 𝑒𝑟𝑟𝑜𝑟𝑒𝑠

𝜆= =
𝑙𝑖𝑏𝑟𝑜 𝑝𝑎𝑔𝑖𝑛𝑎
0.56. DISTRIBUCIONES DISCRETAS GENERALES cxv

El objetivo es calcular 𝑃 (𝑋 = 0), para obtener esta probabilidad de usa la

siguiente instrucción.

dpois(x=0, lambda=0.2)

## [1] 0.8187308
Así 𝑃 (𝑋 = 0) = 0.8187.

0.56. Distribuciones discretas generales

En la práctica nos podemos encontramos con variables aleatorias discretas que
no se ajustan a una de las distribuciones mostradas anteriormente, en esos casos,
es posible manejar ese tipo de variables por medio de unas funciones básicas de
R como se muestra en el siguiente ejemplo.

Ejemplo
El cangrejo de herradura hembra se caracteriza porque su caparazón se adhieren
los machos de la misma especie, en la Figura 6 se muestra una fotografía de este
cangrejo. Los investigadores están interesado en determinar cual es el patrón
de variación del número de machos sobre cada hembra, para esto, se recolectó
una muestra de hembras a las cuales se les observó el color, la condición de la
espina, el peso en kilogramos, el ancho del caparazón en centímetros y el número
de satélites o machos sobre el caparazón, la base de datos está disponible en el
siguiente enlace17 .
1) Encontrar la distribución de probabilidad para la variable Sa que corres-
ponde al número de machos sobre el caparazón de cada hembra.
Primero se debe leer la base de datos usando la url suministrada y luego se
construye la tabla de frecuencia relativa y se almacena en el objeto t1.

url <- '[Link]

crab <- [Link](file=url, header=T)

t1 <- [Link](table(crab$Sa))
t1

##
## 0 1 2 3 4 5
## 0.358381503 0.092485549 0.052023121 0.109826590 0.109826590 0.086705202
## 6 7 8 9 10 11
## 0.075144509 0.023121387 0.034682081 0.017341040 0.017341040 0.005780347
17 [Link]
cxvi DISTRIBUCIONES DISCRETAS

Figura 6: Fotografía del cangrejo de herradura, tomada de

[Link]
tagging

## 12 14 15
## 0.005780347 0.005780347 0.005780347
La anterior tabla de frecuencias relativas se puede representar gráficamente
usando el siguiente código.

plot(t1, las=1, lwd=5, xlab='Número de satélites',

ylab='Proporción')

2) Sea 𝑋 la variable número de satélites por hembra, construir la función

𝐹 (𝑥).
Para construir 𝐹 (𝑥) se utiliza la función ecdf o empirical cumulative density
function, a esta función le debe ingresar el vector con la información de la varia-
ble cuantitativa, a continuación del código usado. En la Figura 8 se muestra la
función de distribución acumulada para para el número de satélites por hembra.
0.56. DISTRIBUCIONES DISCRETAS GENERALES cxvii

0.35

0.30

0.25
Proporción

0.20

0.15

0.10

0.05

0.00

0 1 2 3 4 5 6 7 8 9 10 11 12 14 15

Número de satélites

Figura 7: Función de masa de probabilidad para el número de satélites por

hembra.

F <- ecdf(crab$Sa)
plot(F, las=1, main='')

3) Calcular 𝑃 (𝑋 ≤ 9).
Para obtener esta probabilidad se usa el objeto F que es en realidad una función,
a continuación la instrucción usada.

F(9)

## [1] 0.9595376
Así 𝑃 (𝑋 ≤ 9) = 0.9595.
4) Calcular 𝑃 (𝑋 > 4).
Para obtener esta probabilidad se usa el hecho de que 𝑃 (𝑋 > 4) = 1−𝑃 (𝑋 ≤ 4),
así la instrucción a usar es.

1 - F(4)
cxviii DISTRIBUCIONES DISCRETAS

1.0

0.8

0.6
Fn(x)

0.4

0.2

0.0

0 5 10 15

Figura 8: Función de distribución acumulada para el número de satélites por

hembra.

## [1] 0.2774566
Por lo tanto 𝑃 (𝑋 > 4) = 0.2775.
5) Suponga que el grupo 1 está formado por las hembras cuyo ancho de
caparazón es menor o igual al ancho mediano, el grupo 2 está formado por
las demás hembras. ¿Será 𝐹 (𝑥) diferente para los dos grupos?
Para realizar esto vamos a particionar el vector Sa en los dos grupos de acuerdo
a la nueva variable grupo creada como se muestra a continuacion.

grupo <- ifelse(crab$Wt <= median(crab$Wt), 'Grupo 1', 'Grupo 2')

x <- split(x=crab$Sa, f=grupo)

El objeto x es una lista y para acceder a los vectores allí almacenados usamos
dos corchetes [[]], uno dentro del otro. Luego para calcular 𝐹 (𝑥) para los dos
grupos se procede así:

F1 <- ecdf(x[[1]])
F2 <- ecdf(x[[2]])

Para obtener las dos 𝐹 (𝑥) en la misma figura se usa el código siguiente.
0.56. DISTRIBUCIONES DISCRETAS GENERALES cxix

plot(F1, col='blue', main='', las=1)

plot(F2, col='red', add=T)
legend('bottomright', legend=c('Grupo 1', 'Grupo 2'),
col=c('blue', 'red'), lwd=1)

1.0

0.8

0.6
Fn(x)

0.4

0.2

Grupo 1
0.0 Grupo 2

0 5 10 15

Figura 9: Función de distribución acumulada para el número de satélites por

hembra diferenciando por grupo.

En la Figura 9 se muestran las dos 𝐹 (𝑥), en color azul para el grupo 1 y en

color rojo para el grupo 2. Se observa claramente que las curvas son diferentes
antes de 𝑥 = 9. El hecho de que la curva azul esté por encima de la roja para
valores menores de 9, es decir, 𝐹1 (𝑥) ≥ 𝐹2 (𝑥), indica que las hembras del grupo
1 tienden a tener menos satélites que las del grupo 2, esto es coherente ya que
las del grupo 2 son más grandes en su caparazón.
cxx DISTRIBUCIONES DISCRETAS
Distribuciones continuas

En este capítulo se mostrarán las funciones de R para distribuciones continuas.

0.57. Funciones disponibles para distribuciones

continuas
Para cada distribución continua se tienen 4 funciones, a continuación el listado
de las funciones y su utilidad.

dxxx(x, ...) # Función de densidad de probabilidad, f(x)

pxxx(q, ...) # Función de distribución acumulada hasta q, F(x)
qxxx(p, ...) # Cuantil para el cual P(X <= q) = p
rxxx(n, ...) # Generador de números aleatorios.

En el lugar de las letras xxx se de debe colocar el nombre de la distribución en

R, a continuación el listado de nombres disponibles para las 11 distribuciones
continuas básicas.

beta # Beta
cauchy # Cauchy
chisq # Chi-cuadrada
exp # Exponencial
f # F
gamma # Gama
lnorm # log-normal
norm # normal
t # t-student
unif # Uniforme
weibull # Weibull

Combinando las funciones y los nombres se tiene un total de 44 funciones, por

ejemplo, para obtener la función de densidad de probabilidad 𝑓(𝑥) de una normal

cxxi
cxxii DISTRIBUCIONES CONTINUAS

se usa la función dnorm( ) y para obtener la función acumulada 𝐹 (𝑥) de una

Beta se usa la función pbeta( ).

Ejemplo beta
Considere que una variable aleatoria 𝑋 se distribuye beta con parámetros 𝑎 = 2
y 𝑏 = 5.
1) Dibuje la densidad de la distribución.
La función dbeta sirve para obtener la altura de la curva de una distribución
beta y combinándola con la función curve se puede dibujar la densidad solici-
tada. En la Figura 10 se presenta la densidad, observe que para la combinación
de parámetros 𝑎 = 2 y 𝑏 = 5 la distribución es sesgada a la derecha.

curve(dbeta(x, shape1=2, shape2=5), lwd=3, las=1,

ylab='Densidad')

2.5

2.0
Densidad

1.5

1.0

0.5

0.0

0.0 0.2 0.4 0.6 0.8 1.0

Figura 10: Función de densidad para una 𝐵𝑒𝑡𝑎(2, 5).

2) Calcular 𝑃 (0.3 ≤ 𝑋 ≤ 0.7).

Para obtener la probabilidad o área bajo la densidad se puede usar la función
integrate, los límites de la integral se ingresan por medio de los parámetros
lower y upper. Si la función a integrar tiene parámetros adicionales como en
0.57. FUNCIONES DISPONIBLES PARA DISTRIBUCIONES CONTINUAScxxiii

este caso, éstos parámetros se ingresan luego de los límites de la integral. A

continuación el código necesario para obtener la probabiliad solicitada.

integrate(f=dbeta, lower=0.3, upper=0.7,

shape1=2, shape2=5)

## 0.40924 with absolute error < 4.5e-15

Otra forma de obtener la probabilidad solicitada es restando de 𝐹 (𝑥𝑚𝑎𝑥 ) la
probabilidad 𝐹 (𝑥𝑚𝑖𝑛 ). Las probabilidades acumuladas hasta un valor dado se
obtienen con la función pbeta, a continuación el código necesario.

pbeta(q=0.7, shape1=2, shape2=5) - pbeta(q=0.3, shape1=2, shape2=5)

## [1] 0.40924
De ambas formas se obtiene que 𝑃 (0.3 ≤ 𝑋 ≤ 0.7) = 0.4092.

Recuerde que para distribuciones continuas

𝑃 (𝑎 < 𝑋 < 𝑏) = 𝑃 (𝑎 ≤ 𝑋 < 𝑏) = 𝑃 (𝑎 < 𝑋 ≤ 𝑏) = 𝑃 (𝑎 ≤ 𝑋 ≤ 𝑏)

Ejemplo normal estándar

Suponga que la variable aleatoria 𝑍 se distribuye normal estándar, es decir,
𝑍 ∼ 𝑁 (0, 1).
1) Calcular 𝑃 (𝑍 < 1.45).
Para calcular la probabilidad acumulada hasta un punto dado se usa la función
pnorm y se evalúa en el cuantil indicado, a continuación el código usado.

pnorm(q=1.45)

## [1] 0.9264707
En la Figura 11 se muestra el área sombreada correspondiente a 𝑃 (𝑍 < 1.45).
2) Calcular 𝑃 (𝑍 > −0.37).
Para calcular la probabilidad solicitada se usa nuevamente la función pnorm
evaluada en el cuantil dado. Como el evento de interés es 𝑍 > −0.37, la pro-
babilidad solicitada se obtiene como 1 - pnorm(q=-0.37), esto debido a que
por defecto las probabilidades entregadas por la función pxxx son siempre a
izquierda. A continuación el código usado.
cxxiv DISTRIBUCIONES CONTINUAS

1 - pnorm(q=-0.37)

## [1] 0.6443088

En la Figura 11 se muestra el área sombreada correspondiente a 𝑃 (𝑍 > −0.37).

Otra forma para obtener la probabilidad solicitada sin hacer la resta es usar el
parámetro [Link] para indicar que interesa la probabilidad a la derecha
del cuantil dado, a continuación un código alternativo para obtener la misma
probabilidad.

pnorm(q=-0.37, [Link]=FALSE)

## [1] 0.6443088

3) Calcular 𝑃 (−1.56 < 𝑍 < 2.58).

Para calcular la probabilidad solicitada se obtiene la probabilidad acumulada

hasta 2.58 y de ella se resta lo acumulado hasta -1.56, a continuación el código
usado.

pnorm(q=2.58) - pnorm(-1.56)

## [1] 0.93568

En la Figura 11 se muestra el área sombreada correspondiente a 𝑃 (−1.56 <

𝑍 < 2.58).

4) Calcular el cuantil 𝑞 para el cual se cumple que 𝑃 (𝑍 < 𝑞) = 0.95.

Para calcular el cuantil en el cual se cumple que 𝑃 (𝑍 < 𝑞) = 0.95 se usa la

función qnorm, a continuación el código usado.

qnorm(p=0.95)

## [1] 1.644854

En la Figura 11 se muestra el área sombreada correspondiente a 𝑃 (𝑍 < 𝑞) =

0.95.

El parámetro [Link] es muy útil para indicar si estamos trabajando

una cola a izquierda o una cola a derecha.
0.57. FUNCIONES DISPONIBLES PARA DISTRIBUCIONES CONTINUAScxxv

P(Z < 1.45) P(Z > −0.37)

0.4

0.4
0.3

0.3
Density

Density
0.2

0.2
0.1

0.1
0.0

0.0
−4 −2 0 2 4 −4 −2 0 2 4

x x

P(−1.56 < Z < 2.58) P(Z<q)=0.95

0.4

0.4
0.3

0.3
Density

Density
0.2

0.2
0.1

0.1
0.0

0.0

−4 −2 0 2 4 −4 −2 0 2 4

x x

Figura 11: Área sombreada para los ejemplos.

Ejemplo normal general

Considere un proceso de elaboración de tornillos en una empresa y suponga que
el diámetro de los tornillos sigue una distribución normal con media de 10 𝑚𝑚
y varianza de 4 𝑚𝑚2 .
1) Un tornillo se considera que cumple las especificaciones si su diámetro está
entre 9 y 11 mm. ¿Qué porcentaje de los tornillos cumplen las especifica-
ciones?
Como se solicita probabilidad se debe usar pnorm indicando que la media es
𝜇 = 10 y la desviación de la distribución es 𝜎 = 2. A continuación el código
cxxvi DISTRIBUCIONES CONTINUAS

usado.

pnorm(q=11, mean=10, sd=2) - pnorm(q=9, mean=10, sd=2)

## [1] 0.3829249
2) Un tornillo con un diámetro mayor a 11 mm se puede reprocesar y recu-
perar. ¿Cuál es el porcentaje de reprocesos en la empresa?
Como se solicita una probabilidad a derecha se usa [Link]=FALSE dentro
de la función pnorm. A continuación el código usado.

pnorm(q=11, mean=10, sd=2, [Link]=FALSE)

## [1] 0.3085375
3) El 5 % de los tornillos más delgados no se pueden reprocesar y por lo tanto
son desperdicio. ¿Qué diámetro debe tener un tornillo para ser clasificado
como desperdicio?
Aquí interesa encontrar el cuantil tal que 𝑃 (𝐷𝑖𝑎𝑚𝑒𝑡𝑟𝑜 < 𝑞) = 0.05, por lo tanto
se usa la función qnorm. A continuación el código usado.

qnorm(p=0.05, mean=10, sd=2)

## [1] 6.710293
4) El 10 % de los tornillos más gruesos son considerados como sobredimensio-
nados. ¿cuál es el diámetro mínimo de un tornillo para que sea considerado
como sobredimensionado?
Aquí interesa encontrar el cuantil tal que 𝑃 (𝐷𝑖𝑎𝑚𝑒𝑡𝑟𝑜 > 𝑞) = 0.10, por lo tanto
se usa la función qnorm pero incluyendo [Link]=FALSE por ser una cola a
derecha. A continuación el código usado.

qnorm(p=0.10, mean=10, sd=2, [Link]=FALSE)

## [1] 12.5631
En la Figura 12 se muestran las áreas sombreadas para cada de las anteriores
preguntas.

0.58. Distribuciones continuas generales

En la práctica nos podemos encontramos con variables aleatorias continuas que
no se ajustan a una de las distribuciones mostradas anteriormente, en esos casos,
0.58. DISTRIBUCIONES CONTINUAS GENERALES cxxvii

P(9 < Diámetro < 11) P(Diámetro > 11)

0.20 0.20

0.15 0.15
Densidad

Densidad
0.10 0.10

0.05 0.05

0.00 0.00

4 6 8 10 12 14 16 4 6 8 10 12 14 16

Diámetro Diámetro

P(Diámetro < q) = 5% P(Diámetro > q) = 10%

0.20 0.20

0.15 0.15
Densidad

Densidad

0.10 0.10

0.05 0.05

0.00 0.00

4 6 8 10 12 14 16 4 6 8 10 12 14 16

Diámetro Diámetro

Figura 12: Área sombreada para el ejemplo de los tornillos.

es posible manejar ese tipo de variables por medio de unas funciones básicas de
R como se muestra en el siguiente ejemplo.

Ejemplo
En este ejemplo se retomará la base de datos crab sobre el cangrejo de herradura
hembra presentado en el capítulo anterior. La base de datos crab contiene las
siguientes variables: el color del caparazón, la condición de la espina, el peso
en kilogramos, el ancho del caparazón en centímetros y el número de satélites
o machos sobre el caparazón, la base de datos está disponible en el siguiente
cxxviii DISTRIBUCIONES CONTINUAS

enlace18 .
1) Sea 𝑋 la variable peso del cangrejo, dibuje la densidad para 𝑋.
Para obtener la densidad muestral de un vector cuantitativo se usa la función
density, y para dibujar la densidad se usa la función plot aplicada a un ob-
jeto obtenido con density, a continuación el código necesario para dibujar la
densidad.

url <- '[Link]

crab <- [Link](file=url, header=T)

plot(density(crab$W), main='', lwd=5, las=1,

xlab='Peso (Kg)', ylab='Densidad')

0.15
Densidad

0.10

0.05

0.00

20 25 30 35

Peso (Kg)

Figura 13: Función de densidad 𝑓(𝑥) para el peso de los cangrejos.

En la Figura 13 se muestra la densidad para la variable peso de los cangrejos,

esta densidad es bastante simétrica y el intervalo de mayor densidad está entre
22 y 30 kilogramos.
2) Dibujar 𝐹 (𝑥) para el peso del cangrejo.
Para dibujar la función 𝐹 (𝑥) se usa la función ecdf y se almacena el resultado
18 [Link]
0.58. DISTRIBUCIONES CONTINUAS GENERALES cxxix

en el objeto F, luego se dibuja la función deseada usando plot. A continuación

el código utilizado. En la Figura 14 se presenta el dibujo para 𝐹 (𝑥).

F <- ecdf(crab$W)
plot(F, main='', xlab='Peso (Kg)', ylab='F(x)', cex=0.5, las=1)

1.0

0.8

0.6
F(x)

0.4

0.2

0.0

20 25 30 35

Peso (Kg)

Figura 14: Función acumulada 𝐹 (𝑥) para el peso de los cangrejos.

3) Calcular la probabilidad de que un cangrejo hembra tenga un peso inferior

o igual a 28 kilogramos.
Para obtener 𝑃 (𝑋 ≤ 28) se evalua en la función 𝐹 (𝑥) el cuantil 28 así.

F(28)

## [1] 0.7919075
Por lo tanto 𝑃 (𝑋 ≤ 28) = 0.7919.
4) Dibujar la función de densidad para el peso de los cangrejos hembra dife-
renciando por el color del caparazón.
Como son 4 los colores de los caparazones se deben construir 4 funciones de
densidad. Usando la función split se puede partir el vector de peso de los
cangrejos según su color. Luego se construyen las cuatro densidades usando la
función density aplicada a cada uno de los pesos, a continuación el código.
cxxx DISTRIBUCIONES CONTINUAS

pesos <- split(x=crab$W, f=crab$C)

f1 <- density(pesos[[1]])
f2 <- density(pesos[[2]])
f3 <- density(pesos[[3]])
f4 <- density(pesos[[4]])

Luego de tener las densidades muestrales se procede a dibujar la primera den-

sidad con plot, luego se usa la funció lines para agregar a la densidad inicial
las restantes densidades. En la Figura 15 se muestran las 4 densidades, una por
cada color de caparazón.

plot(f1, main='', las=1, lwd=4,

xlim=c(18, 34),
xlab='Peso (Kg)', ylab='Densidad')
lines(f2, lwd=4, col='red')
lines(f3, lwd=4, col='blue')
lines(f4, lwd=4, col='orange')
legend('topright', lwd=4, bty='n',
col=c('black', 'red', 'blue', 'orange'),
legend=c('Color 1', 'Color 2', 'Color 3', 'Color 4'))

Otra forma para dibujar las densidades es usar el paquete ggplot2 (Wickham
et al., 2021). En la Figura 16 se muestra el resultado obtenido de correr el
siguiente código.

require(ggplot2) # Recuerde que primero debe instalarlo

crab$Color <- [Link](crab$C) # Para convertir en factor

ggplot(crab, aes(x=W)) +
geom_density(aes(group=Color, fill=Color), alpha=0.3) +
xlim(18, 34) + xlab("Peso (Kg)") + ylab("Densidad")

Para aprender más sobre el paquete ggplot2 se recomienda consultar este enla-
ce19 .

19 [Link]
0.58. DISTRIBUCIONES CONTINUAS GENERALES cxxxi

Color 1
0.20 Color 2
Color 3
0.15 Color 4
Densidad

0.10

0.05

0.00

20 25 30

Peso (Kg)

Figura 15: Función de densidad 𝑓(𝑥) para el peso del cangrejo diferenciando
por el color.
cxxxii DISTRIBUCIONES CONTINUAS

0.20

0.15 Color
Densidad

1
2
0.10 3
4

0.05

0.00

20 25 30
Peso (Kg)

Figura 16: Función de densidad 𝑓(𝑥) para el peso del cangrejo diferenciando
por el color y usando ggplot2.
Verosimilitud

En este capítulo se mostrará como usar R para obtener la función de log-

verosimilitud y estimadores por el método de máxima verosimilitud.

0.59. Función de verosimilitud

El concepto de verosimilitud fue propuesto por Fisher (1922) en el contexto de
estimación de parámetros. En la Figura 17 se muestra una fotografía de Ronald
Aylmer Fisher.

Figura 17: Fotografía de Ronald Aylmer Fisher (1890-1962).

La función de verosimilitud para un vector de parámetros 𝜃 dada una muestra

aleatoria 𝑥 = (𝑥1 , … , 𝑥𝑛 )⊤ con una distribución asumida se define usualmente
como:

𝑛
𝐿(𝜃|𝑥) = ∏ 𝑓(𝑥𝑖 |𝜃), (1)
𝑖=1

cxxxiii
cxxxiv VEROSIMILITUD

donde 𝑥𝑖 representa uno de los elementos de la muestra aleatoria y 𝑓(⋅) es la

función de masa/densidad de la distribución de la cual se obtuvo 𝑥.

0.60. Función de log-verosimilitud

La función de log-verosimilitud 𝑙(𝜃|𝑥) se define como el logaritmo de la función
de verosimilitud 𝐿(𝜃|𝑥), es decir
𝑛
𝑙(𝜃|𝑥) = log 𝐿(𝜃|𝑥) = ∑ log 𝑓(𝑥𝑖 |𝜃) (2)
𝑖=1

0.61. Método de máxima verosimilitud para es-

timar parámetros
El método de máxima verosimilitud se usa para estimar los parámetros de una
distribución. El objetivo de este método es encontrar los valores de 𝜃 que ma-
ximizan a 𝐿(𝜃|𝑥) o a 𝑙(𝜃|𝑥), los valores encontrados se representan usualmente
por 𝜃.̂

Asumiendo un modelo estadístico parametrizado por una cantidad fija y

desconocida 𝜃, la verosimilitud 𝐿(𝜃) es la probabilidad de los datos obser-
vados 𝑥 como una función de 𝜃 (Pawitan, 2013). Si la variable de interés
es discreta se usa la probabilidad y si es continua se usa la densidad para
obtener la verosimilitud.

Ejemplo
En este ejemplo vamos a considerar la distribución binomial cuya función de
masa de probabilidad está dada por:

𝑛
𝑓(𝑥) = 𝑃 (𝑋 = 𝑥) = ( )𝑝𝑥 (1 − 𝑝)𝑛−𝑥 , 0 < 𝑝 < 1, 𝑛 ≤ 1, 2, … , 0≤𝑥≤𝑛
𝑥

La distribución binomial anterior tiene sólo un parámetro 𝑝, por lo tanto en este

caso se 𝜃 = 𝑝.
Suponga que se tiene el vector rta que corresponde a una muestra aleatoria de
una distribución binomial con parámetro 𝑛 = 5 conocido.

rta <- c(2, 2, 1, 1, 1, 1, 0, 2, 1, 2,

1, 0, 1, 2, 1, 0, 0, 2, 2, 1)

1) Calcular el valor de log-verosimilitud 𝑙(𝜃) si asumiendo que 𝑝 = 0.30 en la

distribución binomial.
0.61. MÉTODO DE MÁXIMA VEROSIMILITUD PARA ESTIMAR PARÁMETROScxxxv

Para obtener el valor de 𝑙(𝜃) en el punto 𝑝 = 0.30 se aplica la definición dada

en la expresión (2). Como el problema trata de una binomial se usa entonces la
función de masa dbinom evaluada en la muestra rta, el parámetro size como
es conocido se reemplaza por el valor de cinco y en el parámetro prob se cambia
por 0.3. Como interesa la función de log-verosimilitud se debe incluir log=TRUE.
A continuación el código necesario.

sum(dbinom(x=rta, size=5, prob=0.3, log=TRUE))

## [1] -24.55231
Por lo tanto 𝑙(𝜃) = −24.55
2) Construir una función llamada ll a la cual le ingrese valores del parámetro
𝑝 de la binomial y que la función entregue el valor de log-verosimilitud.
La función solicitada tiene un cuerpo igual al usado en el numeral anterior, a
continuación el código necesario para crearla.

ll <- function(prob) sum(dbinom(x=rta, size=5, prob=prob, log=T))

Vamos a probar la función en dos valores arbitrarios 𝑝 = 0.15 y 𝑝 = 0.80 que

pertenezcan al dominio del parámetro 𝑝 de la distribución binomial.

ll(prob=0.15) # Individual para p=0.15

## [1] -25.54468

ll(prob=0.80) # Individual para p=0.80

## [1] -98.45598
El valor de log-verosimilitud para 𝑝 = 0.15 fue de -25.54 mientras que para
𝑝 = 0.80 fue de -98.46.
Vamos ahora a chequear si la función ll está vectorizada y para esto usamos el
código mostrado a continuación y deberíamos obtener un vector con los valores
c(-25.54, -98.56).

ll(prob=c(0.15, 0.80))

## [1] -57.31899
No obtuvimos el resultado esperado, eso significa que nuestra función no está
vectorizada. Ese problema lo podemos solucionar así:
cxxxvi VEROSIMILITUD

ll <- Vectorize(ll)
ll(prob=c(0.15, 0.80))

## [1] -25.54468 -98.45598

Vemos que ahora que cuando se ingresa un vector a la función ll se obtiene un
vector.

Necesitamos que la función ll esté vectorizada para poder dibujarla y

para poder optimizarla.

3) Dibujar la curva log-verosimilitud 𝑙(𝜃), en el eje X debe estar el parámetro

𝑝 del cual depende la función de log-verosimilitud.
En la Figura 18 se presenta la curva solicitada.

curve(ll, lwd=4, col='dodgerblue3',

xlab='Probabilidad de éxito (p)', las=1,
ylab=expression(paste("Probabilidad de éxito (p=", theta, ")"))
)
grid()

4) Observando la Figura 18, ¿cuál esl el valor de 𝑝 que maximiza la función

de log-verosimilitud?
Al observar la Figura 18 se nota que el valor de 𝑝 que maximiza la función
log-verosimilitud está muy cerca de 0.2.
5) ¿Cuál es el valor exacto de 𝑝 que maximiza la función log-verosimilitud?
En R existe la función optimize que sirve para encontrar el valor que minimiza
una función uniparamétrica en un intervalo dado, sin embargo, aquí interesa es
maximimizar la función de log-verosimilitud, por esa razón se construye la fun-
ción minusll que es el negativo de la función ll para así poder usar optimize.
A continuación el código usado.

minusll <- function(x) -ll(x)

optimize(f=minusll, interval=c(0, 1))

## $minimum
## [1] 0.229993
##
## $objective
## [1] 23.3246
0.61. MÉTODO DE MÁXIMA VEROSIMILITUD PARA ESTIMAR PARÁMETROScxxxvii

−50
Probabilidad de éxito (p=θ)

−100

−150

−200

−250

−300

0.0 0.2 0.4 0.6 0.8 1.0

Probabilidad de éxito (p)

Figura 18: Función de log-verosimilitud para el ejemplo sobre binomial.

Del resultado anterior se observa que cuando 𝑝 = 0.23 el valor máximo de log-
verosimilitud es -23.32 (negativo de minusll).

Ejemplo
Suponga que la estatura de una población se puede asumir como una normal
𝑁 (170, 25). Suponga también que se genera una muestra aleatoria de 50 obser-
vaciones de la población con el objetivo de recuperar los valores de la media y
varianza poblacionales a partir de la muestra aleatoria.
La muestra se va a generar con la función rnorm pero antes se fijará una semilla
con la intención de que el lector pueda replicar el ejemplo y obtener la misma
muestra aleatoria aquí generada, el código para hacerlo es el siguiente.

[Link](1235) # La semilla es 1235

y <- rnorm(n=50, mean=170, sd=5)
cxxxviii VEROSIMILITUD

y[1:7] # Para ver los primeros siete valores generados

## [1] 166.5101 163.5757 174.9498 170.5589 170.5710 178.4910 170.2392

1) Construya la función de log-verosimilitud para los parámetros de la normal
dada la muestra aleatoria y.
Abajo se muestra la forma de construir la función de log-verosimilitud.

ll <- function(param) {
media <- param[1] # param es el vector de parámetros
desvi <- param[2]
sum(dnorm(x=y, mean=media, sd=desvi, log=TRUE))
}

Siempre que el interés sea encontrar los valores que maximizan una función
de log-verosimilitud, los parámetros de la distribución deben ingresar a
la función ll como un vector. Esto se debe hacer para poder usar las
funciones de búsqueda optim y nlminb.

2) Dibujar la función de log-verosimilitud.

En la Figura 19 se muestra el gráfico de niveles para la superficie de log-
verosimilitud. De esta figura se nota claramente que los valores que maximizan
la superficie están alrededor de 𝜇 = 170 y 𝜎 = 5.

ll1 <- function(a, b) sum(dnorm(x=y, mean=a, sd=b, log=TRUE))

ll1 <- Vectorize(ll1)
xx <- seq(from=160, to=180, by=0.5)
yy <- seq(from=3, to=7, by=0.5)
zz <- outer(X=xx, Y=yy, ll1)
[Link](x=xx, y=yy, z=zz, nlevels=20,
xlab=expression(mu), ylab=expression(sigma),
color = [Link])

3) Obtenga los valores de 𝜇 y 𝜎 que maximizan la función de log-

verosimilitud.
Para obtener los valores solicitados vamos a usar la función nlminb que es un
optimizador. A la función nlminb se le debe indicar por medio del parámetro
objective la función que queremos optimizar (minimizar); el parámetro start
es un vector con los valores iniciales para comenzar la búsqueda de 𝜇 y 𝜎; los
parámetros lower y upper sirven para delimitar el espacio de búsqueda. A
continuación se muestra el código usado para obtener los valores que minimizan
a minusll, es decir, los valores que maximizan la función de log-verosimilitud.
0.61. MÉTODO DE MÁXIMA VEROSIMILITUD PARA ESTIMAR PARÁMETROScxxxix

7
−150

−200
6
−250

−300
5
σ

−350

4 −400

−450
3
160 165 170 175 180

Figura 19: Gráfico de niveles para la función de log-verosimilitud para el ejem-

plo sobre normal.

minusll <- function(x) -ll(x)

nlminb(objective=minusll, start=c(163, 3.4),
lower=c(160, 3), upper=c(180, 7))

## $par
## [1] 170.338374 5.423529
##
## $objective
## [1] 155.4842
##
## $convergence
## [1] 0
##
## $iterations
## [1] 13
##
## $evaluations
## function gradient
## 16 35
##
## $message
## [1] "relative convergence (4)"
cxl VEROSIMILITUD

De la salida anterior podemos observar que los valores óptimos de 𝜇 y 𝜎 son

170.338 y 5.424 respectivamente, resultado que coincide con lo observado en la
Figura 19 y con los valores reales de simulación de la muestra. Esto indica que el
procedimiento de estimación de parámetros por máxima verosimilitud entrega
valores insesgados de los parámetros a estimar.
Un resultado interesante de la salida anterior es que se reporta el valor mínimo
que alcanza la función minusll, este valor fue de 155.5, por lo tanto, se puede
afirmar que el valor máximo de log-verosimilitud es -155.5.
Otros resultados importantes de la salida anterior son el valor de convergence=0
que indica que la búsqueda fue exitosa; iterations=13 indica que se realizaron
13 pasos desde el punto inicial start hasta las coordenadas de optimización.

En R se tienen dos funciones básicas para optimizar funciones, es decir,

para encontrar los valores que minimizan una función dada. Esas dos fun-
ciones son nliminb y optim. Para optimizar en una sola dimensión se usa
la función optimize.

4) ¿Hay alguna función para obtener directamente el valor que maximiza la

función log-verosimilitud?
La respuesta es si. Si la distribución estudiada es una de las distribuciones
básicas se puede usar la función fitdistr del paquete básico MASS. Esta
función requiere de los datos que se ingresan por medio del parámetro x, y de la
distribución de los datos que se ingresa por medio del parámetro densfun. La
función fitdistr admite 15 distribuciones diferentes para hacer la búsqueda
de los parámetros que caracterizan una distribución, se sugiere consultar la
ayuda de la función fitdistr escribiendo en la consola help(fitdistr). A
continuación el código usado.

require(MASS) # El paquete ya está instalado, solo se debe cargar

res <- fitdistr(x=y, densfun='normal')
res

## mean sd
## 170.3383794 5.4235271
## ( 0.7670026) ( 0.5423527)
El objeto res contiene los resultados de usar fitdistr. En la primer línea están
los valores de los parámetros que maximizan la función de log-verosimilitud, en
la parte de abajo, dentro de paréntesis, están los errores estándar o desviaciones
de éstos estimadores.
Al objeto res es de la clase fitdistr y por lo tanto se le puede aplicar la función
genérica logLik para obtener el valor de la log-verosimilitud. Se sugiere consul-
tar la ayuda de la función logLik escribiendo en la consola help(logLik). A
0.61. MÉTODO DE MÁXIMA VEROSIMILITUD PARA ESTIMAR PARÁMETROScxli

continuación el código para usar logLik sobre el objeto res.

logLik(res)

## 'log Lik.' -155.4842 (df=2)

De esta última salida se observa que el valor coincide con el obtenido cuando se
usó nlminb.

Ejemplo
Generar 𝑛 = 100 observaciones de una gamma con parámetro de forma igual a
2, parámetro de tasa igual a 0.5 y luego responder las preguntas.
1) ¿Cómo se puede generar la muestra aleatoria solicitada?
Para generar la muestra aleatoria (ma) solicitada se fijó la semilla con el objetivo
de que el lector pueda obtener los mismos resultados de este ejemplo.

n <- 100
[Link](12345)
ma <- rgamma(n=n, shape=2, rate=0.5)

2) Asumiendo que la muestra aleatoria proviene de una normal (lo cual es

incorrecto), estime los parámetros de la distribución normal.

fit1 <- fitdistr(x=ma, densfun='normal')

fit1

## mean sd
## 4.3082767 2.8084910
## (0.2808491) (0.1985903)
3) Asumiendo que la muestra aleatoria proviene de una gamma estime los
parámetros de la distribución gamma.

fit2 <- fitdistr(x=ma, densfun='gamma')

fit2

## shape rate
## 2.23978235 0.51987909
## (0.29620136) (0.07702892)
En la salida anterior están los valores estimados de los parámetros de la distri-
bución por el método de máxima verosimilitud, observe la cercanía de éstos con
los verdaderos valores de 2 y 0.5 para forma y tasa respectivamente.
cxlii VEROSIMILITUD

4) Dibuje dos qqplot, uno asumiendo distribución normal y el otro distribu-

ción gamma. ¿Cuál distribución se ajusta mejor a los datos simulados?
Para dibujar el qqplot se usa la función genérica qqplot, recomendamos con-
sultar Hernández (2018) para los detalles de cómo usar esta función. Al usar
qqplot para obtener el qqplot normal y gamma es necesario indicar los valores 𝜃 ̂
obtenidos en el numeral anterior, por eso es que en el código mostrado a continua-
ción aparece mean=4.3083, sd=2.8085 en el qqplot normal y shape=2.23978,
rate=0.51988 en el qqplot gamma.

par(mfrow=c(1, 2))

qqplot(y=ma, pch=19,
x=qnorm(ppoints(n), mean=4.3083, sd=2.8085),
main='Normal Q-Q Plot',
xlab='Theoretical Quantiles',
ylab='Sample Quantiles')

qqplot(y=ma, pch=19,
x=qgamma(ppoints(n), shape=2.23978, rate=0.51988),
main='Gamma Q-Q Plot',
xlab='Theoretical Quantiles',
ylab='Sample Quantiles')

Normal Q−Q Plot Gamma Q−Q Plot

12

12
10

10
Sample Quantiles

Sample Quantiles
8

8
6

6
4

4
2

2
0

0 5 10 0 5 10 15

Theoretical Quantiles Theoretical Quantiles

Figura 20: Gráfico cuantil cuantil normal y gamma para la muestra simulada.

En la Figura 20 se muestran los qqplot solicitados. Se observa claramente que

0.61. MÉTODO DE MÁXIMA VEROSIMILITUD PARA ESTIMAR PARÁMETROScxliii

al asumir normalidad (lo cual es incorrecto), los puntos del qqplot no están
alineados, mientras que al asumir distribución gamma (lo cual es correcto), los
puntos si están alineados. De esta figura se concluye que la muestra ma puede
provenir de una 𝐺𝑎𝑚𝑚𝑎(2.23978, 0.51988).

Para obtener el gráfico cuantil cuantil bajo normalidad se puede usar di-
rectamente la función qqnorm, consultar Hernández (2018) para mayores
detalles.

En este ejemplo se eligió la mejor distribución entre dos candidatas usando

una herramienta gráfica, lo que se recomienda usar algún método menos
subjetivo (cuantitativo) para tomar decisiones.

5) Para comparar modelos se puede utilizar el Akaike information criterion

(𝐴𝐼𝐶) propuesto por Akaike (1974) que sirve para medir la calidad rela-
tiva de los modelos estadísticos, la expresión para calcular el indicador es
𝐴𝐼𝐶 = −2 𝑙 ̂+ 2 𝑑𝑓, donde 𝑙 ̂ corresponde al valor de log-verosimilitud y 𝑑𝑓
corresponde al número de parámetros estimados del modelo. Siempre el
modelo elegido es aquel modelo con el menor valor de 𝐴𝐼𝐶. Calcular el
𝐴𝐼𝐶 para los modelos asumidos normal y gamma.

-2 * logLik(fit1) + 2 * 2 # AIC para modelo normal

## 'log Lik.' 494.3172 (df=2)

-2 * logLik(fit2) + 2 * 2 # AIC para modelo gamma

## 'log Lik.' 466.0479 (df=2)

De los resultados anteriores se concluye que entre los dos modelos, el mejor es
el gamma porque su 𝐴𝐼𝐶 = 466 es el menor de toos los 𝐴𝐼𝐶.

Modelos anidados pueden ser comparados por medio del global deviance
(𝐺𝐷) dado por 𝐺𝐷 = −2 𝑙 ̂ y modelos no anidados por medio del Genera-
lized Akaike information criterion (𝐺𝐴𝐼𝐶) propuesto por Akaike (1983)
y dado por 𝐺𝐴𝐼𝐶 = −2 𝑙 ̂ + ♯ 𝑑𝑓 siendo ♯ el valor de penalidad por cada
parámetro adicional en el modelo; cuando ♯ = 2, el 𝐺𝐴𝐼𝐶 coincide con
el 𝐴𝐼𝐶 y el Schwarz Bayesian criterion (𝑆𝐵𝐶) propuesto por Schwarz
(1978) se dá cuando el valor de penalidad es ♯ = log(𝑛) donde 𝑛 es el
número de observaciones del modelo; siempre el modelo elegido es aquel
modelo con el menor valor de cualquiera de los criterios de información
anteriores.
cxliv VEROSIMILITUD

0.62. Score e Información de Fisher

En esta sección se explican los conceptos y utilidad de la función Score y la
Información de Fisher.

0.62.1. Score e Información de Fisher en el caso univaria-

do.
La función Score denotada por 𝑆(𝜃) se define como la primera derivada de la
función de log-verosimilitud así:

𝜕
𝑆(𝜃) ≡ 𝑙(𝜃)
𝜕𝜃

y el estimador de máxima verosimilitud 𝜃 ̂ se encuentra solucionando la igualdad

𝑆(𝜃) = 0

En el valor máximo 𝜃 ̂ la curva 𝑙(𝜃) es cóncava hacia abajo y por lo tanto la

segunda derivada es negativa, así la curvatura 𝐼(𝜃) se define como

𝜕2
𝐼(𝜃) ≡ − 𝑙(𝜃)
𝜕𝜃2

Una curvatura grande 𝐼(𝜃)̂ está asociada con un gran pico en la función de
log-verosimilitud y eso significa una menor incertidumbre sobre el parámetro 𝜃
(Pawitan, 2013). En particular la varianza del estimador de máxima verosimili-
tud está dada por

𝑉 𝑎𝑟(𝜃)̂ = 𝐼 −1 (𝜃)̂

Ejemplo
Suponga que se desea estudiar una variable que tiene distribución Poisson con
parámetro 𝜆 desconocido. Suponga además que se tienen dos situciones:
Un solo valor 5 para estimar 𝜆.
Cuatro valores 5, 10, 6 y 15 para estimar 𝜆.
Dibujar la función 𝑙(𝜆) para ambos casos e identificar la curvatura.
A continuación el código para evaluar la función 𝑙(𝜆) para cada caso.

# Caso 1
w <- c(5)
ll1 <- function(lambda) sum(dpois(x=w, lambda=lambda, log=T))
ll1 <- Vectorize(ll1)
0.63. MÉTODO DE MÁXIMA VEROSIMILITUD PARA ESTIMAR PARÁMETROS EN MODELOS DE REGRESI

# Caso 2
y <- c(5, 10, 6, 15)
ll2 <- function(lambda) sum(dpois(x=y, lambda=lambda, log=T))
ll2 <- Vectorize(ll2)

En la Figura 21 se muestran las dos curvas 𝑙(𝜆) para cada uno de los casos. De
la figura se observa claramente que cuando se tienen 4 observaciones la curva
es más puntiaguda y por lo tanto menor incertibumbre sobre el parámetro 𝜆 a
estimar.

0
Con 5
Con 5, 10, 6 y 15
−5
l(λ)

−10

−15

−20

0 5 10 15 20 25

Figura 21: Curvas de log-verosimilitud para los dos casos.

0.63. Método de máxima verosimilitud para es-

timar parámetros en modelos de regresión
En esta sección se mostrará como estimar los parámetros de un modelo de
regresión general.
cxlvi VEROSIMILITUD

Ejemplo
Considere el modelo de regresión mostrado abajo. Simule 1000 observaciones del
modelo y use la función optim para estimar los parámetros del modelo.

𝑦𝑖 ∼ 𝑁 (𝜇𝑖 , 𝜎2 ),
𝜇𝑖 = −2 + 3𝑥1 ,
𝜎 = 5,
𝑥1 ∼ 𝑈 (−5, 6).

El código mostrado a continuación permite simular un conjunto de valores con

la estructura anteior.

n <- 1000
x1 <- runif(n=n, min=-5, max=6)
y <- rnorm(n=n, mean=-2 + 3 * x1, sd=5)

El vector de parámetros del modelo anterior es 𝜃 = (𝛽0 , 𝛽1 , 𝜎)⊤ = (−2, 3, 5)⊤ , el

primer elemento corresponde al intercepto, el segundo a la pendiente y el último
a la desviación.

minusll <- function(theta, y, x1) {

media <- theta[1] + theta[2] * x1 # Se define la media
desvi <- theta[3] # Se define la desviación.
- sum(dnorm(x=y, mean=media, sd=desvi, log=TRUE))
}

Ahora vamos a usar la función optim para encontrar los valores que maximizan la
función de log-verosimilitud, el código para hacer eso se muestra a continuación.
En el parámetro par se coloca un vector de posibles valores de 𝜃 para iniciar la
búsqueda, en fn se coloca la función de interés, en lower y upper se colocan
vectores que indican los límites de búsqueda de cada parámetro, los 𝛽𝑘 pueden
variar entre −∞ y ∞ mientras que el parámetro 𝜎 toma valores en el intervalo
(0, ∞). Como la función minusll tiene argumentos adicionales y e x1, estos
pasan a la función optim al final como se muestra en el código.

res1 <- optim(par=c(0, 0, 1), fn=minusll, method='L-BFGS-B',

lower=c(-Inf, -Inf, 0), upper=c(Inf, Inf, Inf),
y=y, x1=x1)

En el objeto res1 está el resultado de la optimización, para explorar los resul-

tados usamos
0.63. MÉTODO DE MÁXIMA VEROSIMILITUD PARA ESTIMAR PARÁMETROS EN MODELOS DE REGRESI

res1

## $par
## [1] -1.904603 3.079600 5.014184
##
## $value
## [1] 3031.209
##
## $counts
## function gradient
## 19 19
##
## $convergence
## [1] 0
##
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
De la salida anterior se observa que el vector de parámetros estimado es 𝛽0̂ =
−1.9046025, 𝛽1̂ = 3.0796 y 𝜎̂ = 5.0141842, se observa también que el valor de
la máxima log-verosimilitud fue de -3031.2092104. Vemos entonces que el vector
estimado está muy cerca del verdadero 𝜃 = (𝛽0 = −2, 𝛽1 = 3, 𝜎 = 5)⊤ .

Cuando se usa optim es necesario decirle que inicie la búsqueda de 𝜃 a

partir de un lugar. Por esa razón se usó par=c(0, 0, 1), esto significa
que la búsqueda inicia en el tripleta 𝛽0 = 0, 𝛽1 = 0 y 𝜎 = 1.

En algunas ocasiones es mejor hacer la búsqueda de los parámetros en el inter-

valo (−∞, ∞) que en una región limitada como por ejemplo (0, ∞) o (−1, 1),
ya que las funciones de búsqueda podrían tener problemas en los bordes de esos
intervalos. Una estrategia usual en este tipo de casos es aplicar una transfor-
mación apropiada al parámetro que tiene el dominio limitado. En el presente
ejemplo 𝜎 sólo puede tomar valores mayores que cero y una transformación de
tipo log podría ser muy útil ya que log relaciona los reales positivos con todos
los reales. La transformación para este problema sería log(𝜎) = 𝛽3 o escrita de
forma inversa 𝜎 = exp(𝛽3 ). El nuevo parámetro 𝛽3 puede variar en (−∞, ∞)
pero al ser transformado por la función exponencial este se volvería un valor
apropiado para 𝜎. Para implementar esta variación lo único que se debe hacer
es modificar la línea 3 de la función minusll como se muestra a continuación:

minusll <- function(theta, y, x1) {

media <- theta[1] + theta[2] * x1
desvi <- exp(theta[3]) # <<<<<---- El cambio fue aquí
- sum(dnorm(x=y, mean=media, sd=desvi, log=TRUE))
cxlviii VEROSIMILITUD

Para hacer la búsqueda se procede de forma similar, abajo el código necesario.

res2 <- optim(par=c(0, 0, 0), fn=minusll, method='L-BFGS-B',

y=y, x1=x1)
res2

## $par
## [1] -1.904609 3.079598 1.612271
##
## $value
## [1] 3031.209
##
## $counts
## function gradient
## 21 21
##
## $convergence
## [1] 0
##
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
De la salida anterior se observa que el vector de parámetros estimado es 𝛽0̂ =
−1.9046094, 𝛽1̂ = 3.0795984 y 𝜎̂ = exp(1.6122706) = 5.0141834, se observa
también que el valor de la máxima log-verosimilitud fue de -3031.2092104. Vemos
entonces que el vector estimado está muy cerca del verdadero 𝜃 = (𝛽0 = −2, 𝛽1 =
3, 𝜎 = 5)⊤ .

EJERCICIOS
1) Al inicio del Capítulo 0.45 se presentó la base de datos sobre medidas
del cuerpo, consulte la explicación sobre la base de datos y responda lo
siguiente.
Si se asume que la edad tiene distribución normal, ¿cuáles son los estima-
dores de máxima verosimilitud para 𝜇 y 𝜎?
Como el histograma para la edad muestra un sesgo a la derecha se podría
pensar que la distribución gamma sería una buena candidata para explicar
las edades observadas. Asumiendo una distribución gamma, ¿cuáles son
los estimadores de máxima verosimilitud para los parámetros?
¿Cuál de los dos modelos es más apropiado para explicar la variable de
interés? Calcule el 𝐴𝐼𝐶 para decidir.
0.63. MÉTODO DE MÁXIMA VEROSIMILITUD PARA ESTIMAR PARÁMETROS EN MODELOS DE REGRESI

2) En el capítulo 0.54 se presentó un ejemplo donde se usó la base de datos

sobre cangrejos hembra. Consulte la explicación sobre la base de datos y
responda lo siguiente.
Suponga que el número de satélites sobre cada hembra es una variable
que se distribuye Poisson. Construya en R la función de log-verosimilitud
𝑙, dibuje la función 𝑙 y encuentre el estimador de máxima verosimilitud de
𝜆.
Repita el ejercicio anterior asumiendo que el número de satélites se distri-
buye binomial negativo.
¿Cuál de los dos modelos es más apropiado para explicar la variable de
interés? Calcule el 𝐴𝐼𝐶 para decidir.
3) Al inicio del Capítulo 0.48 se presentó la base de datos sobre apartamen-
tos usados en Medellín, consulte la explicación sobre la base de datos y
responda lo siguiente.
Dibuje una densidad para la variable área del apartamento.
Describa lo encontrado en esa densidad.
¿Qué distribuciones de 2 parámetros podrían explicar el comportamiento
del área de los apartamentos? Mencione al menos 3.
Para cada una de las distribuciones anteriores dibuje un gráfico de contor-
nos o calor para la función de log-verosimilitud y estime los parámetros
de la distribución elegida.
¿Cuál de los dos modelos es más apropiado para explicar la variable de
interés? Calcule el 𝐴𝐼𝐶 para decidir.
4) Considere el siguiente modelo de regresión.

𝑦𝑖 ∼ 𝐺𝑎𝑚𝑚𝑎(𝑠ℎ𝑎𝑝𝑒𝑖 , 𝑠𝑐𝑎𝑙𝑒𝑖 ),
log(𝑠ℎ𝑎𝑝𝑒𝑖 ) = 3 − 7𝑥1 ,
log(𝑠𝑐𝑎𝑙𝑒𝑖 ) = 3 − 1𝑥2 ,
𝑥1 ∼ 𝑈 (0, 1),
𝑥2 ∼ 𝑃 𝑜𝑖𝑠𝑠𝑜𝑛(𝜆 = 3)

Simule 100 observaciones del modelo anterior.

Escriba el vector de parámetros del problema.
Construya la función minusll para el problema.
Use la función optim para estimar los parámetros del problema.
cl VEROSIMILITUD