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}\)

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

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, 
                 ...)
}

Var(dnorm, sd = 4)

## [1] 16

Var(dunif, min = 10, max = 16)

## [1] 2.999995

Var(dexp, rate = 2)

## [1] 0.25