Duodécima sesión con Rcmdr: funciones de densidad para variables aleatorias continuas.

Definición de funciones en R

En las dos últimas sesiones hemos visto sendos ejemplos de la forma en la que R gestiona una distribución discreta (la binomial, en la 10ª sesión), y una distribución continua (la normal, en la 11ª sesión). Se trata en ambos casos, de distribuciones muy importantes, con nombre propio, por así decirlo, y a las que R conoce bien. Por eso existen funciones específicas como dbinom, pnorm, etc., para trabajar con estas distribuciones en concreto.

Por otra parte, en la novena sesión,  hemos visto como utilizar R para tratar con variables aleatorias discretas genéricas (que no tienen nombre propio, como la binomial), definidas mediante una tabla de densidad de probabilidad.  Y la pregunta es evidente: ¿qué ocurre con las variables aleatorias genéricas de tipo continuo? ¿Cómo las podemos manejar desde R?

Como sabemos, una variable aleatoria continua X se caracteriza mediante una función de densidad, f(x), que permite asignar un valor de probabilidad a cada intervalo de valores de la variable X. Para un intervalo de la forma (a,b), la probabilidad de que X tome valores en ese intervalo se calcula mediante:

\int_a^b f(x)dx

Para poder realizar este tipo de cálculos en R, vamos a aprender dos cosas:

  1. En primer lugar, aprenderemos a definir funciones en R.
  2. A continuación, veremos como calcular sus integrales en un intervalo de valores.

Empecemos con la definición de una función en R. Queremos subrayar desde el principio que, aunque aquí las vamos a usar como funciones de densidad, las funciones se pueden usar para muchas otras cosas en R (veremos ejemplos de esto en próximas entradas).  La definición de funciones (de una o más variables) es una de las características que hacen de R un programa tan flexible y potente.

Supongamos que la función de densidad de nuestra variable aleatoria X es

f(x)=6x(1-x)

en el intervalo (0,1) (más abajo veremos que es realmente una función de densidad). Para definir esta función en R usamos este comando:

f = function(x){ 6*x*(1-x) }

Es decir, f es el nombre que nosotros hemos elegido para la función. A continuación viene el comando function de R, seguido del nombre de la variable (o las variables, puede haber varias separadas por comas). Y finalmente, lo que llamaremos el cuerpo de la función, que es la parte en la que explicamos como calcular el valor de la función a partir de la variable o variables. Aunque no es imprescindible en todos los casos, es bueno acostumbrarse a escribir el cuerpo de la función entre llaves, porque eso hace más fácil la lectura de los ficheros de código que escribimos.

Una vez definida la función, podemos usarla como las otras funciones de R. Por ejemplo, para saber el valor de f(x) en x=1/2  hacemos

f(1/2)

y en este ejemplo se obtiene 1.5 como respuesta. Insistimos, estas funciones que definimos son como las funciones del propio R. Podemos por ejemplo, aplicar la función f(x) a un vector:

datos=c(1,3,-1,2,5)
f(datos)

y se obtiene como respuesta:

> datos=c(1,3,-1,2,5)
> f(datos)
[1]    0  -36  -12  -12 -120

Para dibujar la gráfica de la función de densidad (o de cualquier función de una variable, en realidad) podemos usar el comando curve que vimos en la 11ª sesión. En concreto, para la f(x) que estamos usando de ejemplo, en el intervalo (0,1), hacemos:

curve(f,from=0,to=1,col="red",ylim=range(0, 1.5),lwd=6)

(ver la anterior sesión para la explicación de las opciones) y se obtiene:

Integrales. Cálculo de la media y varianza de una variable aleatoria continua.

Ahora que ya sabemos definir la función de densidad, vamos a pasar a la segunda parte del plan. Aprenderemos a calcular (siempre aproximadamente) la integral de la función en un intervalo. Y como aplicación de esto, vamos a calcular la media y desviación típica de la variable aleatoria definida por esa función de densidad.

La integración se lleva a cabo mediante el comando integrate de R Debemos indicar el nombre de la función, y los extremos del intervalo de integración, mediante las opciones lower y upper. Por ejemplo, la integral de nuestra función de ejemplo f(x) en el intervalo (1/3,2/3) se obtiene mediante:

integrate( f, lower=1/3, upper=2/3)

La respuesta que se obtiene es:

> integrate( f, lower=1/3, upper=2/3 )
0.4814815 with absolute error < 5.3e-15

Como se ve, R no se limita a calcular el valor aproximado de la integral sino que además nos da información sobre la calidad de esa aproximación. Eso está muy bien, pero si queremos utilizar el resultado en otra operación, puede ser un engorro. afortunadamente el remedio es muy  fácil. Sólo hay que añadir una pequeña modificación al final de la llamada al comando integrate:

integrate( f, lower=1/3, upper=2/3 )$value

y ahora se obtiene como salida un número (el valor aproximado de la integral),

> integrate( f, lower=1/3, upper=2/3 )$value
[1] 0.4814815

que podemos usar en otras operaciones. Como segundo ejemplo, vamos a comprobar algo que hemos dejado pendiente. Veamos que la función f(x) que estamos usando es, realmente, una función de densidad en el intervalo (0,1), comprobando que su integral en ese intervalo es 1. Ejecutamos este comando:

integrate(f, lower=0, upper=1)$value

y la respuesta es, en efecto, 1.

¿Y para calcular la media \mu de la variable X cuya densidad es f(x)? Recordemos que, si la densidad de X en (a,b), es una función f(x), entonces la media de X es:

\mu_X=\int_a^b x f(x) dx

Para aplicar esto a la función f(x), como primer paso creamos una función auxiliar, que llamaremos  aux_f (elegimos ese nombre para recordar su papel; podría ser cualquier otro), cuyo valor es el de f(x) multiplicado por x.

aux_f = function(x){ x * f(x) }

Es interesante observar que la definición de esta función incluye una llamada a la propia función f(x), en lugar de copiar directamente su definición. De esa manera, esta función auxiliar sigue siendo válida incluso aunque la función f(x) original cambie,

Y ahora, para calcular la media \mu de f basta con calcular la integral de esta función auxiliar entre 0 y 1.

mu=integrate(aux_f, lower=0, upper=1)$value
mu

Se obtiene \mu=2 como valor de la media.

El siguiente paso lógico es obtener la varianza \sigma (y desviación típica) de X.  Ahora, recordando que la fórmula para la varianza es

\sigma^2_X=\int_a^b (x-\mu)^2 f(x) dx

ya debería estar claro que el primer paso es definir una nueva función auxiliar:

aux2_f = function(x){ (x-mu)^2 * f(x) }

e integrarla entre 0 y 1.

varianza=integrate(aux2_f, lower=0, upper=1)$value
varianza

El valor que se obtiene de la varianza es 0.05=\frac{1}{20}. La desviación típica es simplemente la raíz cuadrada de este número (calculada, por ejemplo, con la función sqrt de R).

Función de distribución (acumulada)

La función de distribución (acumulada) de una variable aleatoria continua X, definida en el intervalo (a,b) por la función de densidad f(x) es

F(x)=P(X\leq x)=\int_a^x f(t)dt

Para aclararlo un poco: si, por ejemplo, (a,b)=(1,5), y queremos calcular F(3), entonces tenemos que integrar la densidad f(x) desde 1 (el principio del intervalo) hasta 3 (el valor en el que calculamos F).

Esta función de distribución F se puede definir en R  de forma sencilla, recurriendo al comando integrate. En lugar de F, que ya se usa para el valor lógico FALSE en R, vamos a llamarla Df en R. Eso sí, primero debemos almacenar en una variable el comienzo del intervalo donde está definida X:

a=0 #extremo inferior del intervalo
Df = function(x){ integrate(f, lower=a, upper=x)$value }

Obsérvese el uso de $value, para obtener el valor de la integral. Ahora podemos pedirle a R que calcule cualquier valor de la función de distribución (acumulada) Df, Por ejemplo, evidentemente Df1)=1 en este ejemplo; pero también se obtiene Df(1/3)=0.2592593, lo cual significa que:

P(X<1/3)=0.2592593

Recordemos que en el caso de la distribución normal, la función de distribución acumulada  (con el mismo sentido que estamos usando aquí) se obtiene mediante pnorm. Y en general, para cualquier distribución “con nombre propio”, el prefijo p indica que calculamos la función de distribución. Lo que hemos aprendido a hacer es calcular esto para nuestras propias variables aleatorias continuas.

Muestreo de variables aleatorias continuas

Cuando, en la novena sesión,  vimos como utilizar R para tratar con variables aleatorias discretas genéricas, terminamos aquella sesión explicando como obtener muestras aleatorias de ese tipo de variables usando la función sample de R.

Desafortunadamente, sample no puede usarse en el caso continuo. Hay, no obstante una forma de obtener una muestra de una variable aleatoria continua dada mediante su función de densidad. El problema es que, para entender porque este método funciona, es necesario conocer más Estadística y más sobre el propio R de lo que hasta ahora hemos visto en el blog. Así que nos vamos a limitar a dar un código que funciona, y a dejar las explicaciones para más adelante.  Una advertencia, eso sí. El cálculo de muestras grandes, con funciones de densidad complicadas, puede llevar bastante tiempo (minutos, horas; en algunos caso puede resultar inviable).

El fichero R que cierra esta entrada contiene, resumidas, la mayor parte de las instrucciones que hemos comentado en esta entrada, incluido el muestreo.

Gracias por la atención.

##############################################################
#
# Estadística, Grado en Biología y Biologia Sanitaria
# Universidad de Alcala
#
# Fichero de instrucciones R para trabajar
# con variables aleatorias continuas genéricas.
##############################################################

#Definción de la función de densidad
f = function(x){ 6*x*(1-x) }

#Intervalo (a,b) de definición de la variable
a=0
b=1

#Comprobación de coherencia
integrate( f, lower=a, upper=b)$value
#El resultado debe ser 1, de lo contrario tenemos un problema.

#Cálculo de la media
aux_f = function(x){ x * f(x) }
mu=integrate(aux_f, lower=a, upper=b)$value
mu

#Cálculo de la varianza
aux2_f = function(x){ (x-mu)^2 * f(x) }
varianza=integrate(aux2_f, lower=a, upper=b)$value
varianza

# y desviación típica
desvTipica=sqrt(varianza)
desvTipica

#Definimos la función de distribución acumulada.
Df = function(x){ integrate(f, lower=a, upper=x)$value }

#Construimos una muestra aleatoria de la variable.

#Número de valores en la muestra
nValores = 1000;
puntosUniformes = runif(nValores);

#El código siguiente genera un vector muestra_f que contiene la muestra aleatoria.
sampling_f = function(x, u) {
return(Df(x) - u);
}

muestra_f = c();
for (i in 1:nValores) {
# find the root within (a,b)
r = uniroot(sampling_f, c(a,b), tol = 0.0001, u = puntosUniformes[i])$root;
muestra_f = c(muestra_f, r);
}
muestra_f

#Puede ser interesante terminar comparando el diagrama de frecuencias de la muestra con la función de densidad original.
h = hist(muestra_f, plot=F)
ylim = range(0, h$density, 0)
hist(muestra_f, freq=F, density=20,ylim=ylim)
curve(f,add=T,col="blue",lwd=4)
Anuncios

Un pensamiento en “Duodécima sesión con Rcmdr: funciones de densidad para variables aleatorias continuas.

  1. Pingback: QQ-plots y funciones de distribución empírica en R. | PostData

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión /  Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión /  Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión /  Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión /  Cambiar )

Conectando a %s