¡Descarga simulacion del software R y más Apuntes en PDF de Estadística Aplicada solo en Docsity!
M. Jesús García-Ligero
Patricia Román Román
SIMULACIÓN CON R
Instalación de R
Página principal de R (si se pone en Google “R”, es la página asociada con “The
R Project for Statistical Computing”):
www. r -project.org/
En la columna izquierda seleccionamos
Download, Packages
CRAN
A continuación buscamos, dentro de CRAN Mirrors, el pais España
Spain
http://cran.es.r-project.org/
Spanish National Research Network,
Madrid
y accedemos a dicha página web.
En la sección
Download and Install R
Precompiled binary distributions of the base system and contributed packages,
Windows and Mac users most likely want one of these versions of R:
Linux
MacOS X
Windows
seleccionamos el sistema operativo de nuestro ordenador.
A continuación, seleccionamos el subdirectorio base
base Binaries for base distribution (managed by Duncan Murdoch)
y descargamos la versión disponible:
Download R 2.10.1 for Windows (32 megabytes)
Todo lo anterior (si se selecciona el sistema operativo Windows) se puede hacer accediendo
directamente a la página
http://cran.es.r-project.org/bin/windows/base
Editor Tinn-R
Se puede descargar por ejemplo de la página http://www.sciviews.org/Tinn-R/
Instalar paquetes
Cargar paquetes
Comenzando con R
Al ejecutar el programa en la interfaz de RGui aparece el símbolo > esperando la
entrada de instrucciones.
El menú principal contiene las pestañas típicas de otras aplicaciones: Archivo, Editar,
Visualizar, Ventanas y Ayuda , junto con una específicas de R: Misc y Paquetes.
Ayuda en R
Para solicitar ayuda sobre un tema general podemos escribir, por ejemplo
help(Uniform)
y obtenemos una ventana de ayuda sobre la distribución uniforme. También podemos
acceder a la misma información a través del menú Ayuda/Funciones R(texto).
Si no se conoce la grafía de una expresión podemos escribir
apropos(“vector”)
y aparecerán las expresiones que contienen el término introducido.
Con la instrucción
help(“vector”)
aparece una ventana de ayuda con información sobre el comando vector.
Operaciones aritméticas
Directamente se pueden realizar operaciones aritméticas después del símbolo >:
+ (suma), - (diferencia), * (producto), / (cociente), ^ (potencia):
[1] 54
[1] 25
[1] -
Si el vector ocupara más de una línea, cada nueva línea empezaría con [n], indicando n
el lugar que ocupa dentro del vector la coordenada que sigue.
Vamos a usar n:m para incluir en un vector los valores n, n+1, n+2, …., m, y con ello
mostraremos lo anteriormente comentado
x<-1: x [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 [19] 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 [37] 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 [55] 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 [73] 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
[91] 91 92 93 94 95 96 97 98 99 100
Se puede cambiar la amplitud de lo que aparece en pantalla con la siguiente función
options(width=40) 1: [1] 1 2 3 4 5 6 7 8 [9] 9 10 11 12 13 14 15 16 [17] 17 18 19 20 21 22 23 24 [25] 25 26 27 28 29 30 31 32 [33] 33 34 35 36 37 38 39 40 [41] 41 42 43 44 45 46 47 48 [49] 49 50 51 52 53 54 55 56
[57] 57 58 59 60 61 62 63 64 [65] 65 66 67 68 69 70 71 72 [73] 73 74 75 76 77 78 79 80 [81] 81 82 83 84 85 86 87 88 [89] 89 90 91 92 93 94 95 96 [97] 97 98 99 100
options(width=60) 1: [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 [14] 14 15 16 17 18 19 20 21 22 23 24 25 26 [27] 27 28 29 30 31 32 33 34 35 36 37 38 39 [40] 40 41 42 43 44 45 46 47 48 49 50 51 52 [53] 53 54 55 56 57 58 59 60 61 62 63 64 65 [66] 66 67 68 69 70 71 72 73 74 75 76 77 78 [79] 79 80 81 82 83 84 85 86 87 88 89 90 91 [92] 92 93 94 95 96 97 98 99 100
Para generar los valores de un vector también se pueden usar las funciones seq() y rep()
Con la instrucción seq(a,b,by=r) o, simplemente seq(a,b,r) , se genera una lista de
números que empieza en a y termina en b , de la forma a, a+r, a+2r, …
seq(1,20) [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 seq(1,20, by=2) [1] 1 3 5 7 9 11 13 15 17 19 seq(20,1) [1] 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 seq(20,1, by=-2)
[1] 20 18 16 14 12 10 8 6 4 2
Con la instrucción seq(a,b,length=r) se generan r números entre a y b , igualmente
espaciados.
seq(4,10, length=8) [1] 4.000000 4.857143 5.714286 6.571429 7.428571 8.285714 9. [8] 10.
Con la instrucción rep(x,r) se genera una lista de r valores todos iguales a x. Se puede
encadenar con la instrucción seq(), aplicarla a vectores (en tal caso, si el número de
repeticiones de cada elemento del vector es distinto, se debe incluir un vector con los
valores de tales repeticiones) y usar la opción each.
rep(3,12) [1] 3 3 3 3 3 3 3 3 3 3 3 3 rep(seq(2,20,by=2),2) [1] 2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16 18 20 rep(c(1,4),c(3,2)) [1] 1 1 1 4 4 rep(c(1,4), each=3) [1] 1 1 1 4 4 4
Vectores
Ya se ha visto como asignar valores a un vector
x<-c(3,6,8,9)
Par ver el contenido del vector basta escribir su nombre
x [1] 3 6 8 9
También se ha visto como usar el símbolo : se puede usar para crear secuencias
crecientes (o decrecientes) de valores. Por ejemplo
Numerosde5a20<-5: Numerosde5a [1] 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Numerosde20a5<-20: Numerosde20a [1] 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5
Se pueden concatenar vectores:
x<-c(1,2) y<-c(3,4) z<-c(x,y) x [1] 1 2 y
Suma de vectores
y<-c(6,7,8,9,10) x+y [1] 7 9 11 13 15
Notemos que las operaciones con vectores se hacen elemento a elemento
Potencia de un vector (eleva al cuadrado cada elemento del vector)
x^ [1] 1 4 9 16 25
Raiz cuadrada (hace la raíz cuadrada de cada elemento del vector)
sqrt(x) [1] 1.000000 1.414214 1.732051 2.000000 2.
Producto de dos vectores (elemento a elemento)
¡Cuidado, con el producto!
x<-c(1,2,3,4,5) y<-c(6,7,8,9,10) x*y [1] 6 14 24 36 50
¡Y la potencia! (Eleva cada elemento del primer vector al correspondiente del segundo)
x^y [1] 1 128 6561 262144 9765625
Cuando los vectores tienen dimensiones diferentes, el de menor dimensión se extiende
repitiendo los valores desde el principio, aunque da un mensaje de aviso
x<-c(1,2,3,4,5) y<-c(1,2,3,4,5,6) x+y [1] 2 4 6 8 10 7 Mensajes de aviso perdidos In x + y : longitud de objeto mayor no es múltiplo de la longitud de uno menor x^y [1] 1 4 27 256 3125 1 Mensajes de aviso perdidos In x^y : longitud de objeto mayor no es múltiplo de la longitud de uno menor
Funciones con vectores
Asignar nombres a los elementos de un vector con la función names()
x<-c(3.141592, 2.718281) names(x)<-c("pi","e")
x pi e 3.141592 2. x["pi"] pi
La función length() devuelve la dimensión de un vector
x<-c(1,2,3,4,5) length(x) [1] 5
Las funciones sum() y cumsum() proporcionan la suma y sumas acumuladas de las
componentes de un vector
sum(x) [1] 15 cumsum(x) [1] 1 3 6 10 15
Las funciones max() y min() proporcionan el valor máximo y mínimo de las
componentes de un vector
max(x) [1] 5 min(x) [1] 1
Las funciones mean(), median(), var() y sd() proporcionan la media, mediana,
cuasivarianza y cuasidesviación típica, respectivamente, de las componentes de un
vector
mean(x) [1] 3 median(x) [1] 3 var(x) [1] 2. sd(x) [1] 1.
Las funciones prod() y cumprod() proporcionan el producto y productos acumulados
de las componentes de un vector
prod(x) [1] 120 cumprod(x) [1] 1 2 6 24 120
La función quantile() proporciona los cuartiles
Variables y operaciones lógicas
Supongamos que tenemos el vector x<-1:20, podemos ver qué elementos son iguales a
un valor con el operador = =
x<-11: x [1] 11 12 13 14 15 16 17 18 19 20 x== [1] FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE
La instrucción anterior devuelve un vector de valores lógicos (o booleanos). FALSE
indica que la condición chequeada (en este caso, ser igual a 15) no se cumple mientras
que TRUE indica que sí se cumple.
De manera análoga se pueden hacer comparaciones con < , > , <= , >= , ¡= (para distinto)
x<-11: x [1] 11 12 13 14 15 16 17 18 19 20 x== [1] FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE x< [1] TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE x> [1] FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE x<= [1] TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE x>= [1] FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE x!= [1] TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE
Si aplicamos la función sum() al vector lógico resultante, R fuerza TRUE al valor
númerico 1 y FALSE al 0, de manera que obtenemos el número de elementos del vector
que cumplen la condición impuesta:
> sum(x==15)
[1] 1
> sum(x<15)
[1] 4
> sum(x>15)
[1] 5
> sum(x<=15)
[1] 5
> sum(x>=15)
[1] 6
> sum(x!=15)
[1] 9
Ejemplo
Vamos a calcular la media y mediana de un vector y el número de valores que están por
debajo de la media y de la mediana
x<-c(1,5,7,9,3,5,6,2,4,7,5,6,9,8,6,2,6,1,4) mean(x) [1] 5. median(x) [1] 5 sum(x<mean(x)) [1] 10 sum(x<median(x)) [1] 7 length(x) [1] 19 x==median(x)
[1] FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE FALSE [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
sum(x==median(x)) [1] 3
Veamos el efecto de colocar un vector lógico entre los corchetes de índice de un vector
z<-1: z[c(TRUE,FALSE,TRUE,FALSE,TRUE)] [1] 1 3 5 z<-3: z[c(TRUE,FALSE,TRUE,FALSE,TRUE)] [1] 3 5 7
Las instrucción x[x<r] devuelve los valores del vector x que verifican la condición
impuesta x<r.
En el siguiente ejemplo se usa la función runif() que genera valores de una distribución
uniforme y que se comentará con posterioridad.
x<-runif(10) x [1] 0.1976741 0.3802387 0.5560528 0.8759437 0.2752229 0.7574220 0. [8] 0.5145545 0.2543729 0. x[x<0.2] [1] 0. x[x<0.4] [1] 0.1976741 0.3802387 0.2752229 0.2543729 0.
Podemos usar operadores lógicos dentro de los corchetes con condiciones compuestas, |
es el operador lógico “o” y & es “y”
x [1] 0.1976741 0.3802387 0.5560528 0.8759437 0.2752229 0.7574220 0. [8] 0.5145545 0.2543729 0. x[(x<0.2)|(x>0.8)] [1] 0.1976741 0. sum((x<0.2)|(x>0.8))
[1] 50.
mean(y) [1] 50.
Nota: Es útil ver los objetos que tenemos guardados en cada momento para no usar
nombres que contengan valores que queramos conservar. Se hace con la función
objects()
Gráficos en R
Diagrama de barras
Sea
x<-c(1,1,1,1,1,2,2,3,3,3,5,6,6,7,7,7)
la función table() tabula los datos en x y da lugar a
table(x) x 1 2 3 5 6 7 5 2 3 1 2 3
La instrucción barplot(table(x)) muestra el diagrama de barras asociado
1 2 3 5 6 7
0
1
2
3
4
5
Histograma
Se realiza con la función hist().
x<-runif(100) hist(x)
Histogram of x
x
Frequency
0.0 0.2 0.4 0.6 0.8 1.
0
5
10
15
Scatterplots
Se realizan con la instrucción plot(x,y)
x<-runif(100) y<-runif(100) plot(x,y)
0.0 0.2 0.4 0.6 0.8 1.
x
y
QQ-plots
Instrucción qqplot(x,y).
En los siguientes ejemplos se utilizan algunas funciones como rnorm() , rt() que
generan valores aleatorios de distribuciones normales y t de Student y que se comentan
con posterioridad.
x<-rnorm(1000) y<-rnorm(1000) qqplot(x,y,main="x e y con la misma distribución normal " )
-4 -2 0 2
0
5
10
15
20
x N(0,1), c la exponencial de una N(0,1)
x
c
Otras instrucciones de R que necesitaremos
El proceso de realización de simulaciones es a menudo muy repetitivo, por lo que
usaremos la sentencia for() que permite repetir una cierta operación un número de
veces.
Su sintaxis es
for ( Nombre_del indice in vector ){comandos}
de forma que al ejecutar dicha instrucción crea una variable llamada Nombre_
del_indice igual a cada uno de los elementos del vector, en secuencia. Para cada valor se
ejecuta todos los comandos incluidos entre llaves, considerándolos como un único
comando. Si sólo hay que ejecutar un solo comando las llaves no son necesarias, pero es
conveniente incluirlas para una mayor claridad.
Cuando se construye un vector a partir de comandos usando for() (o cualquier otra
función de programación) es necesario definir previamente dicho vector. La instrucción
Nombre<-numeric(length=r) o, simplemente, Nombre<-numeric(r ) crea un objeto de
tipo numérico de longitud r.
x<-numeric(10) x [1] 0 0 0 0 0 0 0 0 0 0
Así, si queremos obtener los primeros 12 números de Fibonacci, se procede de la
siguiente forma
Fibonacci_12<-numeric(12) Fibonacci_12[1]<-Fibonacci_12[2]<- for(i in 3:12){Fibonacci_12[i]<-Fibonacci_12[i-2]+Fibonacci_12[i-1]} Fibonacci_ [1] 1 1 2 3 5 8 13 21 34 55 89 144
GENERACIÓN DE NÚMEROS PSEUDOALEATORIOS
Método congruencial multiplicativo
Ejemplo
Generar 50 números pseudoaleatorios a partir del generador congruencial multiplicativo
xn = 171 xn-1 (mod 30269)
un = xn/
con semilla 27218.
Código R
random.number<-numeric(50) random.seed<- for (j in 1:50)
- {random.seed<-(171*random.seed)%%
- random.number[j]<-random.seed/
- }
random.number [1] 0.76385080 0.61848756 0.76137302 0.19478675 0.30853348 0.
[7] 0.82757937 0.51607255 0.24840596 0.47741914 0.63867323 0. [13] 0.44391952 0.91023820 0.65073177 0.27513297 0.04773861 0. [19] 0.92470845 0.12514454 0.39971588 0.35141564 0.09207440 0. [25] 0.34751726 0.42545178 0.75225478 0.63556774 0.68208398 0. [31] 0.81766824 0.82126929 0.43704780 0.73517460 0.71485678 0. [37] 0.12722587 0.75562457 0.21180085 0.21794575 0.26872378 0. [43] 0.75195745 0.58472364 0.98774324 0.90409330 0.59995375 0. [49] 0.24754700 0.
O bien
x<-numeric(50) semilla<- x[1]=(171semilla)%% for(i in 2:50){x[i]=(171x[i-1])%%30269} NumerosAleatorios<-x/ NumerosAleatorios [1] 0.76385080 0.61848756 0.76137302 0.19478675 0.30853348 0. [7] 0.82757937 0.51607255 0.24840596 0.47741914 0.63867323 0. [13] 0.44391952 0.91023820 0.65073177 0.27513297 0.04773861 0.
„qunif‟ proporciona la función de cuantiles
„runif‟ genera valores aleatorios.
Uso:
dunif(x, min=0, max=1, log = FALSE)
punif(q, min=0, max=1, lower.tail = TRUE, log.p = FALSE)
qunif(p, min=0, max=1, lower.tail = TRUE, log.p = FALSE)
runif(n, min=0, max=1)
Argumentos:
x,q: vector de cuantiles.
p: vector de probabilidades.
n: número de observaciones. Si no se especifica se toma igual a 1
min,max: extremos inferior y superior del intervalo que determina la distribución.
Deben ser finitos. Si no se especifican se toman los valores por defecto 0 y 1. Para el
caso min = max = u, el caso degenerado X = u se considera, aunque como no tiene
función de densidad, la función „dunif‟ devuelve „NaN‟ (condición de error).
log, log.p: son valores lógicos; si son TRUE, las probabilidades p se dan como
probabilities log(p).
lower.tail: es un valor lógico; si es TRUE (por defecto), las probabilidades son
P[X ≤x], en otro caso P[X > x].
La forma de uso más habitual para generar números pseudoaleatorios de una
distribución U(a,b) es
runif(n, min=a, max=b)
y para una U(0,1)
runif(n)
Ejemplo
Generar 10 números aleatorios en el intervalo (0,1) y 15 en el intervalo (-1,2)
runif(10) [1] 0.31264669 0.59918637 0.37935813 0.97517055 0.79116964 0. [7] 0.02291407 0.48761985 0.52772211 0. runif(15,min=-1,max=2) [1] 1.96646178 -0.53018219 0.99929876 1.64215831 0.
[7] -0.68915028 1.86967652 -0.04384054 1.32471288 -0.30430071 -
[13] 0.81967883 1.70135362 0.
Si se quiere ejecutar la función anterior, pero partiendo de una semilla concreta (para
garantizar el mismo resultado en cualquier ejecución), se usará la función set.seed()
antes de la anterior
runif(5) [1] 0.3126467 0.5991864 0.3793581 0.9751705 0. runif(5) [1] 0.47277362 0.02291407 0.48761985 0.52772211 0. set.seed(32789) runif(5) [1] 0.3575211 0.3537589 0.2672321 0.9969302 0. set.seed(32789) runif(5) [1] 0.3575211 0.3537589 0.2672321 0.9969302 0.
Ejercicios propuestos
1. Genera 20 números pseudoaleatorios usando
xn = 172 xn-1 (mod 30307)
con semilla inicial x 0 = 17218.
random.number<-numeric(20) random.seed<- for (j in 1:20)
- {random.seed<-(172*random.seed)%%
- random.number[j]<-random.seed/
- }
random.number [1] 0.71656713 0.24954631 0.92196522 0.57801828 0.41914409 0. [7] 0.95882139 0.91727984 0.77213185 0.80667833 0.74867192 0. [13] 0.71019896 0.15422180 0.52614907 0.49764081 0.59421916 0. [19] 0.37954928 0.
Otra forma posible
x0<- x<-numeric(20) x[1]<-(172x0)%% for(i in 2:20){ x[i] <-(x[i-1]172)%%30307} x/ [1] 0.71656713 0.24954631 0.92196522 0.57801828 0.41914409 0.
