Calcular integrales en R

En R podemos utilizar la función integrate para calcular una integral definida, incluso si es una integral impropia. Probemos con una integral muy sencilla:

\[ \int_{0}^{1}xdx = x^2 \bigr\rvert_{0}^{1} = 1^2 - 0^2 = 1 \]

¡Esta será la única integral que resolvamos analíticamente! :)

Comprobemos entonces nuestro cálculo con R utilizando integrate. Veamos los argumentos más importantes de esta función:

  1. f: Es otra función. Es necesario que el primer argumento de esta función sea numérico.
  2. lower: Límite inferior de la integral. Puede ser Inf o -Inf.
  3. upper: Límite superior de la integral. Puede ser Inf o -Inf.

Para integrar funciones matemáticas cuyo argumento sea x, lo más fácil es escribir function(x) seguido de la expresión de la función. Por ejemplo, function(x) x sería el equivalente a \(f(x) = x\) ¡Sencillo! Veámoslo en acción.

integrate(function(x) x, lower = 0, upper = 1)
## 0.5 with absolute error < 5.6e-15

Probemos ahora con una integral impropia con un límite infinito

\[ \int_{0}^{+\infty}e^{-x}dx \]

integrate(function(x) exp(-x), lower = 0, upper = Inf)
## 1 with absolute error < 5.7e-05

Y ahora probemos con una integral impropia con límites finitos pero en la que la función se vaya a infinito en uno de los límites.

\[ \int_{0}^{3}\frac{1}{\sqrt{3-x}}dx \]

integrate(function(x) 1 / sqrt(3 - x), lower = 0, upper = 3)
## 3.464102 with absolute error < 5e-13

Calcular esperanzas

Vamos a crear la función expected_value para calcular la esperanza de cualquier variable aleatoria sólo con tener su función de densidad. Como sabemos, calculamos la esperanza o valor esperado de una variable aleatoria continua con la siguiente integral:

\[ E[X] = \int_{-\infty}^{+\infty}xf_X(x)dx \]

La función integrate devuelve una lista con 5 componentes: value, abs.error, subdivisions, message y call. El elemento que más nos interesa es value, que es el valor de la integral.

expected_value <- function(.density_function, ...) {
  integral <- integrate(function(x) x * .density_function(x, ...), 
                        lower = -Inf, 
                        upper = Inf)
  
  integral$value
}

El argumento .density_function sería la función de densidad, mientras que pasamos el argumento especial ... llamado ellipsis para indicar los parámetros de dicha función de densidad.

Distribución Normal

Decimos que \(X\) sigue una distribución normal de parámetros \(\mu\) y \(\sigma^2\) y lo denotamos como \(X \sim N(\mu, \sigma^2)\) si su función de densidad es la siguiente:

\[ f_X(x) = \frac{1}{\sigma\sqrt{2\pi}}e^{- \frac{1}{2} \left( \frac{x-\mu}{\sigma} \right)^2}I_{\mathbb{R}}(x) \]

Y sabemos que

  • \(E[X] = \mu\)
  • \(Var[X] = \sigma^2\)

En R tenemos por defecto la función dnorm para esta función de densidad, así que no tenemos que escribirla nosotros. Por ejemplo, si \(\mu=10\) y \(\sigma=5\), calculamos el valor esperado de \(X\) como

expected_value(dnorm, mean = 10, sd = 5)
## [1] 10

¡Coincide con lo que nos dice la teoría!

Distribución uniforme continua

Decimos que \(X\) sigue una distribución continua de parámetros \(a\) y \(b\) y lo denotamos como \(X \sim U_{[a, b]}\) si su función de densidad es la siguiente:

\[ f_X(x) = \frac{1}{b - a}I_{[a, b]}(x) \]

Y sabemos que

  • \(E[X] = \frac{a + b}{2}\)
  • \(Var[X] = \frac{(b - a)^2}{12}\)

El intervalo asociado con esta distribución puede ser cerrado, abierto o semi-abierto en cualquiera de los dos lados.

En R tenemos por defecto la función dunif para esta función de densidad, así que tampoco tenemos que escribirla. Por ejemplo, si \(a=10\) y \(b=16\), calculamos el valor esperado de \(X\) como

expected_value(dunif, min = 10, max = 16)
## [1] 12.99998

En este caso, tenemos un pequeño error en el resultado. Esto es el resultado de que R utilice métodos numéricos para calcular la integral. Podríamos ajustar algunos argumentos en integrate para solucionar esto. Sin embargo, en la práctica, una aproximación como esta podría ser suficientemente buena y además está bueno que veas en estos ejemplos que los cálculos de las computadoras pueden tener errores aunque hagas todo bien.

Distribución Exponencial

Decimos que \(X\) sigue una distribución exponencial de parámetro \(\lambda\) y lo denotamos por \(X \sim exp(\lambda)\) si su función de densidad es la siguiente:

\[ f_X(x) = \lambda e^{-\lambda x} I_{[0, +\infty)}(x) \]

Y sabemos que

  • \(E[X] = \frac{1}{\lambda}\)
  • \(Var[X] = \frac{1}{\lambda^2}\)

En R tenemos por defecto la función dexp para esta función de densidad, así que no tenemos que escribirla nosotros. Por ejemplo, si \(\lambda=10\), calculamos el valor esperado de \(X\) como

expected_value(dexp, rate = 10)
## [1] 0.1

Verificar si una función es función de densidad

Podríamos también crear la función is_density para comprobar que .density_function sea efectivamente una función de densidad. Esto lo hacemos verificando que la integral desde \(-\infty\) hasta \(+\infty\) sea igual a 1. Formalmente, decimos que la función \(f_X(x)\) es una función de densidad si y sólo si:

\[ \int_{-\infty}^{+\infty}f_X(x)dx = 1 \]

Vamos a crear nuestra función en R entonces. Como debemos tener en cuenta que la igualdad puede ser aproximada, vamos a usar all.equal para la verificación, ya que esta devolverá TRUE en estos casos, admitiendo por defecto un error cercano a \(1.5 * 10^{-8}\) según la documentación de la función.

Además, como all.equal devuelve un texto en caso de que los valores no sean iguales, vamos a verificar que devuelva TRUE utilizando la función isTRUE.

is_density <- function(.func, ..., .tolerance = 1.5e-8) {
  integral <- integrate(function(x) .func(x, ...), 
                        lower = -Inf, 
                        upper = Inf)
  
  equals_one = isTRUE(all.equal(integral$value, 1, tolerance = .tolerance))
  
  if (!equals_one) {
    cat(sprintf("Warning: function is not a density function. Integral equals %0.9f", 
                integral$value))
  }
  
  equals_one
}

Distribución normal: dnorm

is_density(dnorm, mean = 10)
## [1] TRUE

Distribución uniforme continua: dunif

is_density(dunif, min = 10, max = 16, .tolerance = 1e-5)
## [1] TRUE

Nuevamente, está atento con estos resultados basados en cálculo numérico. Con la tolerancia por defecto de all.equal, nuestra función nos diría que esto no es una función de densidad.

Distribución \(t\)-student con dt

is_density(dt, df = 2)
## [1] TRUE

Una función que no es función de densidad:

is_density(function(x) 2 * dnorm(x))
## Warning: function is not a density function. Integral equals 2.000000000
## [1] FALSE

Incorporar la verificación en el cálculo de la esperanza

Estaría bueno entonces que para calcular la esperanza primero verifiquemos si la función que indicamos efectivamente es una función de densidad.

expected_value <- function(.density_function, ...) {
  if(!is_density(.density_function, ..., .tolerance = 1e-5)) {
    return(NA)
  }
  
  integral <- integrate(function(x) x * .density_function(x, ...), 
                        lower = -Inf, 
                        upper = Inf)
  integral$value
}

expected_value(function(x) 2 * dnorm(x))
## Warning: function is not a density function. Integral equals 2.000000000
## [1] NA

Calcular la esperanza de una función de una variable aleatoria

Otra cosa que podríamos mejorar en expected_value es calcular la esperanza de cualquier función \(g(X)\) de \(X\). Esto serviría para poder obtener la varianza, por ejemplo.

Dada una variable aleatoria continua \(X\) con función de densidad \(f_X(x)\), la esperanza de \(g(X)\) se define de la siguiente manera:

\[ E[g(X)] = \int_{-\infty}^{+\infty}g(x)f_X(x)dx \]

Por defecto, la función de la que queremos obtener la esperanza sería la función identidad, \(g(x) = x\), que en R es identity.

expected_value <- function(.density_function, ..., .g = identity) {
  if(!is_density(.density_function, ..., .tolerance = 1e-5)) {
    return(NA)
  }
  
  integral <- integrate(function(x) .g(x) * .density_function(x, ...), 
                        lower = -Inf, 
                        upper = Inf)
  
  integral$value
}

expected_value(dunif, min = 10, max = 16)
## [1] 12.99998
expected_value(dnorm, .g = function(x) x^2)
## [1] 1
expected_value(dnorm, .g = function(x) (x - expected_value(dnorm, sd = 4))^2, sd = 4)
## [1] 16

Cálculo de la varianza

Utilizando nuestra nueva expected_value, podríamos crear fácilmente la función Var (con la V mayúscula) para calcular la varianza. Recordemos que la varianza de una variable aleatoria es

\[ Var[X] = E[(x - E[X])^2] \]

De este modo, \(Var[X]\) sería un caso especial de la esperanza de una función de \(X\), con \(g(x) = (x-E[X])^2\).

Var <- function(.density_function, ...) {
  mu <- expected_value(.density_function, ...)
  expected_value(.density_function, 
                 .g = function(x) (x - mu)^2, 
                 ...)
}
Distribución normal
Var(dnorm, sd = 4)
## [1] 16
Distribución uniforme continua
Var(dunif, min = 10, max = 16)
## [1] 2.999995
Distribución exponencial
Var(dexp, rate = 2)
## [1] 0.25