1 Introducción

R es un lenguaje gratuito para computación estadística (Linux, MacOS, Windows)
1991 creado por Ross Ihaka y Robert Gentleman, en 1993 se anuncia la pormera versión publica
1995 se convierte en software gratuito.
Es la versión no comercial de S desarollado por John Chambers y otros en Bell Labs (1976), inicialmente porgramdo en Fortran, en 1988 se escribe en C, 1998 sale la versión 4.
Ayuda en la red
Versiones anuales
“Menos amigable que MATLAB pero más que C++ y Fortran”
“I have written software professionally in perhaps a dozen programming languages, and the hardest language for me to learn has been R,” writes consultant John D. Cook

1.1 Instalación

En el orden inidicado

1.2 RStudio

Console window (command window)
Editor Window
History, Environment
Files, Plots, Package
Documentación de RStudio

Ventanas de R

1.3 Directorio de trabajo y ayuda

Direcotrio de trabajo comando

getwd() # obtiene la carpeta en la cual se esta trabajando
setwd("~/Dropbox/UNAM/Econometria")
setwd("C:/Sharon/Documents/RProjects")
dir() #despliega los archivos del directorio

Directorio de trabajo con ambiente gráfico Session-Set Working Directory
Ayuda en console

help.start()
help(matrix) # aplica también para funciones
?matrix # aplica también para funciones 
apropos("matrix") 
example("matrix")
args(rnorm) # conocer los argumentos de una función

Ayuda en la web
- https://rseek.org/
Comentarios

# Ignorado por el programa, por linea

Librerias disponibles

library()

Instalar paquetes, actualizar

update.packages() # Actualizar todos los paquetes
install.packages("mcsm")
update.packages("mcsm")

Librerias disponibles

library()

Cargar libreria disponible

library(MASS)
library(coda)
library(mcsm)

2 Datos simples (clases otomicas)

null (por ejemplo, c())
logicos (TRUE, FALSE, o T, F) la constante NA (Not Available)
numéricos (3.0, 0.132)
enteros
complejos (-1-2i,complex(1,-1,2))
caracteres (“hola”, “mundo”)

con la instrucció attributes() se pueden conocer los atributos de la clase especifica.

2.1 Operaciones directas con datos simples

Se pueden realizar operaciones directamente en la consola

2 ^ 2 - 3.5
5 < 45
help(Arithmetic)
??Relational

Operadores Aritmeticos Operadores de relación

??Relational

Asignación

a = 2 ^ 2 - 3.5 # No estándar
b <- 23;
c <- b + a
c
c * 7
a <- a + 4.5
a
ls() #despliega los nombres de las variables generadas
print(c) #iprime el valor de

Funciones con “datos simples”

rnorm(a)
rnorm(n = a, mean = 3, sd = 5)
x <- rnorm(n = 10000, mean = 3, sd = 5)
plot(x)
hist(x)
x <- runif(n = 10000)
plot(x)
hist(x)

Funciones útiles

help(Trig)
help(log)
help(Special)

Funciones Trigonométricas Funciones Logaritmo y exponencial Funciones especiales

2.2 Algunas instrucciones útiles

Borrar consola

cat("\014")

Borrar variables

rm(a)
rm(list=ls())

Constantes

pi
exp(1)
Inf
1 / 0
NaN

2.3 Scripts

File-New File-R Script
source(‘Script1.R’)
Desde la interfaz RStudio

3 Datos compuestos (clases)

vector
list
matrix
array
factor
time-series
data.frame

3.1 Vectores

c(…) datos del mismo tipo “concatenate”

a <- c(1 + 0i, 2 + 4i) # números complejos
a
a <- 9 : 29 # genera una secuencia de numeros desde el 9 al 29
a
a <- vector("numeric", length = 10) # función vector()
a
a <- c(-1, 5, 45, 9.5)
a
names(a) <- c("nombre1", "nombre2", "nombre3", "nombre4") # podemos nombrar los elementos
a

b <- c(TRUE, FALSE)
b
c <- c() # vector vacio

b <- c("hola", "mundo", "adios", "cruel")
b

Coersión

b <- c(1.7, "hola") # caracteres
b
b <- c(T, 2.3) # numerico
b
class(b) # nos indica que clase de dato es
as.logical(b) # coersión explicita
as.complex(b)
b <- c("hola", TRUE, FALSE, F) # caracteres
b

Elementos

b <- c("comida", "comida", "zeta", "camion", "zeta", "comida")
b[1] # cuenta desde el elemento 1
b[2:4]
b[b < "Tito"]
u <- b < "Tito"
u
b[u]

Operaciones

a
2 * a
a * a
vec2 <- seq(from=0, to=1, by=0.3)
vec2
a + vec2 
vec2 > 0.5
vec2 == 0
vec2 * a
vec2 / a

Funciones aplicadas a vectores.

a
d <- sum(a)
d
t(a)
mean(a)
sqrt(a)
exp(a)

3.1.1 Ejercicio

Cálcula la siguiente expresión para $n\in\{10,100,1000,10000\}$ con $x_i\sim\mathcal{U}(0,1)$

\[S_n=\frac{5}{n\sqrt{2\pi}}\sum_{i=1}^ne^{-\frac{(5x_i-2.5)^2}{2}}\]

n <- 10000
x <- runif(n)
Sn <- (5 / (n * (sqrt(2 * pi)))) * sum(exp(-(1/2) * (5 * x - 2.5) ^ 2))
Sn

## [1] 0.981423

3.2 Listas

Elementos de diferentes tipos

L <- list(uno = 1, dos = c(1,2), tres = seq(0, 1, length = 5), cuatro = c("hola", "adios"))
L
names(L)
L$tres+10
L[2] # es una sublista de L 
L[[2]] # es un vector númerico, semanticamente igual a $
class(L[2])
class(L[[2]])
L["cuatro"]
L[["cuatro"]]
L[c(1,3)]
L[[c(3,3)]]
L[[3]][[3]]

3.3 Matrices

Son vectores con un atributo de dimensión, este atributo es un vector de enteros de longitud 2 (nrow, ncol)

Definición con número de columnas o filas

A <- matrix(data = c(1, 2 , 1, 4, 1, 6), ncol = 3)
A
B <- matrix(1:9, byrow= TRUE, nrow = 3)
B 
x <- matrix(rnorm(15), nrow = 3)
x
m <- matrix(nrow = 2, ncol = 2)
m
dim(m)
attributes(m)
dimnames(m) <- list(c("fil1", "fil2"), c("col1", "col2")) #nombramos a las filas y columnas
m

Definición concatenando vectores o motrices, por filas o columnas

x <- c(1, 3, 5)
y <- c(2, 4, 6)
m1 <- cbind(x, y)
m1
m2 <- rbind(x, y)
m2

Elementos

A
A[2, 3]
A[1,] # como vector
class(A[1,])
A[1, , drop= FALSE] # como matriz
class(A[1, , drop= FALSE])
A[, 2]
A[1:2, 1]
A[[5]]

Operaciones

A <- matrix(1:4, 2, 2); B <- matrix(1:4,2,2)
A / 5
A * B
A / B
A %*% B

Funciones

mean(A)

Otras funciones asociadas
rownames
colnames
rowSums
colSums

3.4 Factors

Son utilizados para representar datos categoricos (nominales: sin orden, ordinales: con orden), pueden ser ordenados o no.

x <-  factor(c("si", "si", "no", "si", "no"))
x # los levels son rodenados alfabeticamente
table(x)
unclass(x)
x <- factor(c("si", "si", "no", "si", "no"), levels = c("si", "no")) # indicar el orden de los levels, incluso podrian ser cosas distintas
x
summary(x)

# Create speed_vector
speed_vector <- c("medium", "slow", "slow", "medium", "fast")

# Convert speed_vector to ordered factor vector
factor_speed_vector <- factor(speed_vector,
                              ordered = TRUE,
                              levels = c("slow", "medium", "fast"))

# Print factor_speed_vector
factor_speed_vector

## [1] medium slow   slow   medium fast  
## Levels: slow < medium < fast

summary(factor_speed_vector)

##   slow medium   fast 
##      2      2      1

factor_speed_vector[2] > factor_speed_vector[5]

## [1] FALSE

3.5 Data frames

Matriz con nombres es las columnas, se pueden pensar tambien como listas donde cada elemento tiene la misma longitud, los nombres (elementos de la lista) van en las columnas y la cantidad de elementos en las filas. A diferecnia de las matrices pueden almancenar diferentes tipos de datos. Las instrucciones read.table() y read.csv() regresan data frames. De ser posible se puede convertir un data frame a una matriz utilizando data.matrix. Tiene 3 atributos: names, row.names, class

df <- data.frame(x = c(11, 12, 14), y = c(19, 20, 21), z = c(10, 9, 7))
df
data.matrix(df)
x <- data.frame(c1 = 1:3, c2 = c("a", "b", "c"))
x
attributes(x)
nrow(x)
ncol(x)

head(df)
tail(df)
str(df)

Elementos

df$z
df[1,2]
df["x"]
df[["x"]]
df[1:2, ]
df[, 2]
df[1:2, "y"]
subset(df, y == 20)

Funciones

mean(df$z)

3.6 Datos no disponibles

NA (datos faltantes) o NaN para operaciones matemáticas indefinidas. Un NaN es un NA pero no el converso no es cierto

j = c(1, 2, NA)
j
max(j, na.rm = TRUE)
j
is.na(j)
is.nan(j)
nodisponible <- is.na(j)
j[nodisponible]
j[!nodisponible]

4 Gráficas

library(datasets)
hist(airquality$Ozone) 

library(datasets)
with(airquality, plot(Wind, Ozone))

library(datasets)
airquality <- transform(airquality, Month = factor(Month))
boxplot(Ozone ~ Month, airquality, xlab = "Month", ylab = "Ozone (ppb)")
with(airquality, plot(Wind, Ozone))
title(main = "Ozone and Wind in New York City") ## Add a title
with(airquality, plot(Wind, Ozone))
title(main = "Ozone and Wind in New York City") ## Add a title
with(airquality, plot(Wind, Ozone, main = "Ozone and Wind in New York City",
 type = "n"))
with(subset(airquality, Month == 5), points(Wind, Ozone, col = "blue"))
with(subset(airquality, Month != 5), points(Wind, Ozone, col = "red"))
legend("topright", pch = 1, col = c("blue", "red"), legend = c("May", "Other Months"))

with(airquality, plot(Wind, Ozone, main = "Ozone and Wind in New York City",
 pch = 20))
model <- lm(Ozone ~ Wind, airquality)
abline(model, lwd = 2)

par(mfrow = c(1, 3), mar = c(4, 4, 2, 1), oma = c(0, 0, 2, 0))
with(airquality, {
 plot(Wind, Ozone, main = "Ozone and Wind")
 plot(Solar.R, Ozone, main = "Ozone and Solar Radiation")
 plot(Temp, Ozone, main = "Ozone and Temperature")
 mtext("Ozone and Weather in New York City", outer = TRUE)
})

?par

with(faithful, plot(eruptions, waiting)) ## Make plot appear on screen device
title(main = "Old Faithful Geyser data") ## Annotate with a title

plot(rnorm(100), type="l", col="gold")

help(plot)
t=rnorm(100)
plot(t, type="l", ylim=range(t), lwd=3, col=rgb(1,0,0,0.3))
lines(t, type="s", lwd=2,col=rgb(0.3,0.4,0.3,0.9))
points(t, pch=20, cex=4, col=rgb(0,0,1,0.3))

4.1 Gráficas con Lattice

Es un sistema independiente al sistema base, por lo que se tiene que cargar la libreria, algunas de sus funciones son

xyplot
bwplot
histogram

A diferecnia del sistema base, las gráficas se construyen con una sola instrucción. Un ejemplo xyplot(y ~ x | f * g, data) aqui y depende de x y f, g son variables condicionales para relaizar la grafica, estas son opcionales, data representa el conjunto de datos donde se gráficaran las variables.

library(datasets)
xyplot(Ozone~Wind, data = airquality)
xyplot(Ozone~Wind, data = airquality,pch = 8,col = "red", main = "Big Apple Data")

4.2 Gráficas con ggplot2

Es una librería que se puede instalar, la documentación se encuentra en Página de ggplo2. Combina los dos sistemas anteriores, base y el lattice, la función

qplot

similar a plot del sistema base, pero sin indicar tanto detalle similar a como funciona lattice, la otra función que es más versatil es ggplot

#necesitamos importar los datos mpg de la libreria ggplot
#qplot(displ, hwy, data=mpg, color = drv)
#qplot(displ, hwy, data=mpg, color = drv, geom = c("point","smooth"))
#qplot(drv, hwy, data = mpg, geom = "boxplot")
#qplot(drv, hwy, data = mpg, geom = "boxplot", color = manufacturer)
#qplot(displ,hwy, data = mpg, facets = .~drv)
#qplot(hwy, data = mpg, facets = drv~.,binwidth=2)

A “grammar” of graphics means that ggplot2 contains building blocks with which you can create your own graphical objects. What are these basic components of ggplot2 plots? There are 7 of them.

Obviously, there’s a DATA FRAME which contains the data you’re trying to plot. Then the AESTHETIC MAPPINGS determine how data are mapped to color, size, etc.The GEOMS (geometric objects) are what you see in the plot (points, lines,shapes) and FACETS are the panels used in conditional plots. You’ve used these or seen them used in the first ggplot2 (qplot) lesson.

There are 3 more. STATS are statistical transformations such as binning,quantiles, and smoothing which ggplot2 applies to the data. SCALES show what coding an aesthetic map uses (for example, male = red, female = blue). Finally, the plots are depicted on a COORDINATE SYSTEM. When you use qplot these were taken care of for you.

# cargar los mpg de la libreria ggplot2
#qplot(displ, hwy, data = mpg, geom=c("point","smooth"), facets = .~drv)

#g <- ggplot(mpg, aes(displ,hwy) )
#g+geom_point()
#g+geom_point()+geom_smooth()
#g+geom_point()+geom_smooth(method = "lm")
#g+geom_point()+geom_smooth(method = "lm")+facet_grid(.~drv)
#g+geom_point()+geom_smooth(method = "lm")+facet_grid(.~drv)+ggtitle("Swirl Rules!")
#g+geom_point(color="pink", size=4, alpha=1/2)
#g+geom_point(size=4, alpha=1/2,aes(color=drv))
#g+geom_point(aes(color= drv))+labs(title="Swirl Rules!")+labs(x="Displacement", y="Hwy Mileage")
#g+geom_point(aes(color=drv), size=2, alpha=1/2)+geom_smooth(size=4,linetype=3,method="lm",se=FALSE)
#g+geom_point(aes(color=drv))+theme_bw(base_family="Times")
#plot(myx, myy, type="l", ylim=c(-3,3))

#aqui no plaica nada diecto son datos que se craron
#g <- ggplot(data = testdat, aes(x=myx,y=myy))
# g+geom_line()
#g+geom_line()+ylim(-3,3)#no es lo mejor
#g+geom_line()+coord_cartesian(ylim=c(-3,3))

#######
#g <- ggplot(data=mpg, aes(x=displ, y=hwy, color=factor(year)))
#g+geom_point()
#g+geom_point()+facet_grid(drv~cyl,margins=TRUE)
#g+geom_point()+facet_grid(drv~cyl,margins=TRUE)+geom_smooth(method="lm", se=FALSE,size=2, color="black")
#g+geom_point()+facet_grid(drv~cyl,margins=TRUE)+geom_smooth(method="lm", se=FALSE,size=2, color="black")+labs(x="Displacement", y="Highway Mileage", title="Swirl Rules!")

#Se tiene que cargar diamonds from ggplot
#qplot(price, data = diamonds)
#qplot(price, data = diamonds, binwidth =18497/30, fill=cut)
#qplot(price, data = diamonds, geom="density")
#qplot(price, data = diamonds, geom="density", color=cut)

##### Ahora scatterplot

#qplot(carat, price, data = diamonds)
#qplot(carat, price, data = diamonds, shape=cut)
#qplot(carat, price, data = diamonds, color=cut)
#qplot(carat, price, data = diamonds, color=cut)+geom_smooth(method = "lm")
#qplot(carat, price, data = diamonds, color=cut, facets = .~cut)+geom_smooth(method = "lm")

###### ahora con ggplot

#g <- ggplot(diamonds, aes(depth, price))
# g+geom_point(alpha=1/3)
#cutpoints <- quantile(diamonds$carat,seq(0,1,length=4),na.rm=TRUE)
#diamonds$car2 <- cut(diamonds$carat,cutpoints)
#g <- ggplot(diamonds, aes(depth, price))
#g+geom_point(alpha=1/3)+facet_grid(cut~car2)
#g+geom_point(alpha=1/3)+facet_grid(cut~car2)+geom_smooth(method = "lm", size=3, color="pink")
#ggplot(diamonds, aes(carat, price))+geom_boxplot()+facet_grid(.~cut)

5 Clustering

5.1 Clustering jerarquico

Por ejemplo en $\mathbb{R}^2$ si se tienen diferentes puntos medimos las distancias entre ellos, y agrupamos en clusters de 2 componentes a los más cercanos, ya que terminamos con esto obsevamos quien esta más cerca hacemos una segunda agrupación, puede que las agupaciones nuevas sea la union de dos anteriores, o la grupación anterior se le agrega un elemento, podemos visualizar estas agrupaciones creando un frame con dos columnas (las coordenads x, y) y utilizar las siguientes instrucciones

‘dist’
‘hclust’
‘plot’

el plot genera un dendrograma donde podemos hacer cortes para observar las agrupaciones de acurdo al valor del corte.

Se pueden medir las distancias entre los clusters, como:

Tomadno la mayor distancia entre pares que pertenenzcan a un mismo cluster
Definiendo un centroide para cada cluster.

Podemos utilizar los heatmap para visualizar los datos con un mapa de colores, en particular estos datos tienen que estar en un matriz

5.2 K-means clustering

Los datos se parten en k-grupos, tal que la suma de los cuadrados de las distancias a los centroides es minimizada. Se requieren tres cosas:

Una metrica
Numro de clusters
Centroides de los clusters

Al finalizar de ordenar de esta manera, regresa los centroides y los elementos asociados a cada centroide.

Procedemos de la siguiente manera, una vez dados el número de clusters y la localización, calculamos la distacia a cada centroide para cada punto, despues nos quedamos con la minima de estas, y así asiganos la pertenencia al cluster. Ahora se recalculan los centroides, reasignandolos a los centroides de los puntos del cluster especifico.

6 Compresión de datos

Se busca el conjunto más pequeño de varibles no correlacionadas, y que expliquen la mayor cantidad de la varianza de los datos.

Análisis de componentes principales (PCA)
Descomposición de valeres singulares (SVD).

\[X= UDV^t\] Las matrices $U$ y $V$ tienen columnas ortogonales (no correlacionadas). Las columnas de de $U$ son lo vectores singualres izquierdos y las columnas de de $V$ son lo vectores singualres derechos de X. $D$ es una matriz diagonal, los valores de dela diagonal en $D$ son los valores singulares de $X$ .

7 Rregresión lineal

library(datasets)
data(mtcars)

model.uno <- lm(mpg ~ wt, mtcars)
confint(model.uno)

8 Leer y escribir datos

8.1 Leer datos

read.table, para leer datos tabulares, con los igueintes argumento
- file, el nombre del archivo o conección
- header, si tiene o no una linea de cabecera
- sep, una cadena que indica como estan separadas als columnas
- colClasses, un vector de caracteres indicando la clase de cada columna de datos
- nrows, el numero de filas en el conjunto de datos
- comment.char, vecto de caracters indicando un comentario
- skip, el numero de lineas a omitir desde desde el principio.
- stringsAsFactors, si las caracteres deben de leerse como factores.
read.csv, para leer datos tabulares separados por comas
readLines, para leer lineas de un archivo de texto.
source, para leer códigos de R (inverso a dump)
dget, para leer códigos de R (inverso a dput)
load, para leer de un worspace guardado
unserialize, Para leer un objeto de R en forma binaria.

if(!file.exists("data")){
dir.create("data")  
}
# Datos tipo CSV
fileUrl <- "https://data.baltimorecity.gov/api/views/dz54-2aru/rows.csv?accesType=DOWNLOAD"
download.file(fileUrl,destfile="./data/cameras.csv", method="curl")

# Datos tipo EXCEL
fileUrl <- "https://data.baltimorecity.gov/api/views/dz54-2aru/rows.xlsx?accesType=DOWNLOAD"
download.file(fileUrl,destfile="./data/cameras.xlsx", method="curl")

cameraData <- read.table("./data/cameras.csv", sep=",", header = TRUE)
head(cameraData)
#cameraData <- read.csv("./data/cameras.csv")
#head(cameraData)

library(xlsx)
#cameraData <- read.xlsx("./data/cameras.xlsx", sheetIndex=1, header = TRUE)
#head(cameraData)

8.2 Escribir datos

write.table
writeLines
dump
dput
save
seriealize

9 Programación

9.1 if-else

a=4

if(a==4){
  print("hola mundo")
}

if(a<5){
d=3 
d
}else{
  d=4
  d
}

if (a<2){
  print(a)
} else if (a==5){
  print(a+3)
} else{
  print("Ño")
}

9.2 for

Las instrucciones break y next funcionan igual que en otros lenguajes de programcación

contador <- 0  
c <- runif(1000)
for (i in c) {

  if(i < 0.5)
  {
    contador <- contador + 1
  }
}
contador

for(i in 1:10) 
  { print(i)}

x<-c("a","b","c")
for(variable in x)
  print(variable)

x <- matrix(1:6, 2, 3)
x
for(i in seq_len(nrow(x))){
  for(j in seq_len(ncol(x))){
    if(j == 2) next
    print(x[i, j])
  }
}

9.3 while

contador=0
numero=1
while(numero<=30){
  if(numero%%5==0){
    contador<-contador+1
  }
  numero<-numero+1
}
contador

9.4 repeat

contador<-0

repeat
  {

    if(contador==2)
  {
    break
  } 
  else
  {
      print(contador)
    contador<- contador+1
  }
  
    }

9.5 Funciones

desv = function(x){
  sqrt(var(x))
}
x<- 1:10 
desv(x) 

fun1 = function(arg1, arg2=1 )
 {
 w = arg1 ^ 2
 return(arg2 + w)
}

fun1(2)
fun1(arg1 = 3, arg2 = 5)

fun1(3,5)


f <- function(a, b) {
 a^2
}
f(2)



# ... se ussan para pasar argumentos variables, los argumentos listados despues de ... deben 
# de ser usados con nombre 

f = function(...)
{
datos = list(...)
medias = lapply(datos,mean)    # lapply aplica una función sobre una lista
varianzas = lapply(datos,var)
maximos = lapply(datos,max)
minimos = lapply(datos,min)

for(i in 1:length(datos))
{
cat("Distribución ",i,": \n")   # La función cat es para visualizar cosas
cat("media: ",medias[[i]],"varianza: ",varianzas[[i]],"maximo: ",maximos[[i]],"minimo: ",minimos[[i]],"\n")
cat("------------------------------------------------\n")
}

}

f(c(1,2),c(1,3,5,7),c(-1,2,-5,6,9))

x = rnorm(100)
y = runif(50)
f(x,y)

make.power <- function(n) {
 pow <- function(x) {
 x^n
 }
 pow
}

cube <- make.power(3)
square <- make.power(2)
cube(3)
square(3)
ls(environment(cube))

zz <- 10
f <- function(x) {
 zz <- 2
 zz^2 + g(x)
}
g <- function(x) {
 x*zz
}
f(3)

9.6 Paquetes

install.packages(“ggvis”) library, require

search() despliega que paquetes están cargado en el ambiente

9.7 lapply sapply

lapply toma una lista como entrada, aplica una función a cada elemento de la lista, regresa una lista de la misma dimensión que la original. la “l” es de lista, lapply(lista,func)
sapply toma una lista como entrada, aplica una función a cada elemento de la lista, regresa una “simplificación” (de ser posible, si no es como lappy) de la aplicación de la función (respecto al tipo de dato que regresa la función) de la misma dimensión que la original.

x <- list(a = 1:5, b = rnorm(10)) 
res <- lapply(x, mean)
str(res)
res <- unlist(res) # se convierte la lista en un vector
str(res) 
funcion1 <- function(a, b) {
  a*b
}
resn <- lapply(x, funcion1, b = 2)
resn
print("este es sapply")
sapply(x, mean)
x<-1:4
lapply(x, runif)
 x <- list(a = matrix(1:4, 2, 2), b = matrix(1:6, 3, 2))
lapply(x, function(elt) elt[,1])

9.8 apply mapply

apply es usado para evaluar una función (usialmente anonima) sobre los margenes de un arreglo. NO es más rápida que un ciclo pero es una sola linea

str(apply)
x <- matrix(rnorm(200), 20, 10)
Sumacolumnas = apply(x, 2, sum)
Sumafilas = apply(x, 1, sum)
Sumacolumnas
Sumafilas
apply(x, 1, quantile, probs = c(0.25, 0.75))

mapply aplica una función recorriendo diferentes valores dados para sus argumentos, sirve para hacer vectorización (en el sentido de generar un vector a partir de una instrucción)

mapply(rep, 1:4, 4:1)  # list(rep(1, 4), rep(2, 3), rep(3, 2), rep(4, 1))

## [[1]]
## [1] 1 1 1 1
## 
## [[2]]
## [1] 2 2 2
## 
## [[3]]
## [1] 3 3
## 
## [[4]]
## [1] 4

9.9 vapply tapply

vapply permite especificar el formato de salida (por loq ue lleva tres argumentos), si el formato de salida especificado no coincide con lo que se ha obtenido entonces
tapply aplica una función a un determinado subconjunto de los datos, diciendo a que indices lo aplicas, por eso lleva tres argumentos como mínimo.

x <- c(rnorm(10), runif(10), rnorm(10, 1))
x
f<-gl(3,10)
f
tapply(x, f, mean)
tapply(x, f, range)

9.10 split

toma objetos y los separa en grupos determinados por un factor o lista uan lista de factores

x <- c(rnorm(10), runif(10), rnorm(10, 1))
x
f<-gl(3,10)
split(x, f)
lapply(split(x, f), mean)

library(datasets)
head(airquality)

##   Ozone Solar.R Wind Temp Month Day
## 1    41     190  7.4   67     5   1
## 2    36     118  8.0   72     5   2
## 3    12     149 12.6   74     5   3
## 4    18     313 11.5   62     5   4
## 5    NA      NA 14.3   56     5   5
## 6    28      NA 14.9   66     5   6

s <- split(airquality, airquality$Month)
lapply(s, function(x) colMeans(x[, c("Ozone", "Solar.R", "Wind")]))

## $`5`
##    Ozone  Solar.R     Wind 
##       NA       NA 11.62258 
## 
## $`6`
##     Ozone   Solar.R      Wind 
##        NA 190.16667  10.26667 
## 
## $`7`
##      Ozone    Solar.R       Wind 
##         NA 216.483871   8.941935 
## 
## $`8`
##    Ozone  Solar.R     Wind 
##       NA       NA 8.793548 
## 
## $`9`
##    Ozone  Solar.R     Wind 
##       NA 167.4333  10.1800

sapply(s, function(x) colMeans(x[, c("Ozone", "Solar.R", "Wind")]))

##                5         6          7        8        9
## Ozone         NA        NA         NA       NA       NA
## Solar.R       NA 190.16667 216.483871       NA 167.4333
## Wind    11.62258  10.26667   8.941935 8.793548  10.1800

sapply(s, function(x) colMeans(x[, c("Ozone", "Solar.R", "Wind")], na.rm = TRUE))

##                 5         6          7          8         9
## Ozone    23.61538  29.44444  59.115385  59.961538  31.44828
## Solar.R 181.29630 190.16667 216.483871 171.857143 167.43333
## Wind     11.62258  10.26667   8.941935   8.793548  10.18000

9.11 Fechas y tiempo

x <- as.Date("1970-01-01")
x
unclass(x)
unclass(as.Date("1970-01-02"))

x <- Sys.time()
x
p <- as.POSIXlt(x)
names(unclass(p))
p$sec

x <- Sys.time()
x
unclass(x)
# x$sec   error
p <- as.POSIXlt(x)
p$sec

datestring <- c("January 10, 2012 10:40", "December 9, 2011 9:10")
x <- strptime(datestring, "%B %d, %Y %H:%M")
x

x <- as.Date("2012-01-01")
y <- strptime("9 Jan 2011 11:34:21", "%d %b %Y %H:%M:%S")
x <- as.POSIXlt(x)
x-y

x <- as.Date("2012-03-01") 
y <- as.Date("2012-02-28")
x-y
x <- as.POSIXct("2012-10-25 01:00:00")
y <- as.POSIXct("2012-10-25 06:00:00", tz = "GMT")
y-x

10 Leer Datos

10.1 Leer Bases de datos con MySQL

instalar MySQl Mac os, descargar datos de la universida de califiornia enlace

install.packages("RMySQL")

library(DBI)

## Warning: package 'DBI' was built under R version 3.4.4

library(RMySQL)

## Warning: package 'RMySQL' was built under R version 3.4.4

ucscDb <- dbConnect(MySQL(), user="genome",
                    host="genome-mysql.cse.ucsc.edu")
result <- dbGetQuery(ucscDb, "show databases;")
dbDisconnect(ucscDb)

## [1] TRUE

result

##               Database
## 1              acaChl1
## 2              ailMel1
## 3              allMis1
## 4              allSin1
## 5              amaVit1
## 6              anaPla1
## 7              ancCey1
## 8              angJap1
## 9              anoCar1
## 10             anoCar2
## 11             anoGam1
## 12             anoGam3
## 13             apaSpi1
## 14             apaVit1
## 15             apiMel1
## 16             apiMel2
## 17             aplCal1
## 18             aptFor1
## 19             aptMan1
## 20             aquChr2
## 21             araMac1
## 22             ascSuu1
## 23             balAcu1
## 24             balPav1
## 25             bisBis1
## 26             bosTau2
## 27             bosTau3
## 28             bosTau4
## 29             bosTau5
## 30             bosTau6
## 31             bosTau7
## 32             bosTau8
## 33           bosTauMd3
## 34             braFlo1
## 35             bruMal2
## 36             bucRhi1
## 37             burXyl1
## 38             caeAng2
## 39             caeJap1
## 40             caeJap4
## 41              caePb1
## 42              caePb2
## 43              caePb3
## 44             caeRem2
## 45             caeRem3
## 46             caeRem4
## 47            caeSp111
## 48             caeSp51
## 49             calAnn1
## 50             calJac1
## 51             calJac3
## 52             calMil1
## 53             canFam1
## 54             canFam2
## 55             canFam3
## 56             capCar1
## 57             carCri1
## 58             cavPor3
## 59                 cb1
## 60                 cb3
## 61                 cb4
## 62                ce10
## 63                ce11
## 64                 ce2
## 65                 ce4
## 66                 ce6
## 67             cerSim1
## 68             chaVoc2
## 69             cheMyd1
## 70             chlSab2
## 71             chlUnd1
## 72             choHof1
## 73             chrPic1
## 74             chrPic2
## 75                 ci1
## 76                 ci2
## 77                 ci3
## 78             colLiv1
## 79             colStr1
## 80             corBra1
## 81             corCor1
## 82             cotJap2
## 83             criGri1
## 84         criGriChoV1
## 85             cucCan1
## 86             danRer1
## 87            danRer10
## 88            danRer11
## 89             danRer2
## 90             danRer3
## 91             danRer4
## 92             danRer5
## 93             danRer6
## 94             danRer7
## 95             dasNov3
## 96             dipOrd1
## 97             dirImm1
## 98                 dm1
## 99                 dm2
## 100                dm3
## 101                dm6
## 102                dp2
## 103                dp3
## 104            droAna1
## 105            droAna2
## 106            droEre1
## 107            droGri1
## 108            droMoj1
## 109            droMoj2
## 110            droPer1
## 111            droSec1
## 112            droSim1
## 113            droSim2
## 114            droVir1
## 115            droVir2
## 116            droYak1
## 117            droYak2
## 118            eboVir3
## 119            echTel1
## 120            echTel2
## 121            egrGar1
## 122            equCab1
## 123            equCab2
## 124            eriEur1
## 125            eriEur2
## 126            eurHel1
## 127            falChe1
## 128            falPer1
## 129            felCat3
## 130            felCat4
## 131            felCat5
## 132            felCat8
## 133            ficAlb2
## 134                fr1
## 135                fr2
## 136                fr3
## 137            fulGla1
## 138            gadMor1
## 139            galGal2
## 140            galGal3
## 141            galGal4
## 142            galGal5
## 143            galGal6
## 144            galVar1
## 145            gasAcu1
## 146            gavSte1
## 147             gbMeta
## 148            geoFor1
## 149                 go
## 150           go080130
## 151           go140213
## 152           go150121
## 153           go180426
## 154            gorGor3
## 155            gorGor4
## 156            gorGor5
## 157            haeCon2
## 158            halAlb1
## 159            halLeu1
## 160            hetBac1
## 161            hetGla1
## 162            hetGla2
## 163               hg16
## 164               hg17
## 165               hg18
## 166               hg19
## 167        hg19Patch10
## 168        hg19Patch13
## 169               hg38
## 170        hg38Patch11
## 171            hgFixed
## 172          hgcentral
## 173 information_schema
## 174            latCha1
## 175            lepDis1
## 176            letCam1
## 177            loaLoa1
## 178            loxAfr3
## 179            macEug1
## 180            macEug2
## 181            macFas5
## 182            manPen1
## 183            melGal1
## 184            melGal5
## 185            melHap1
## 186            melInc2
## 187            melUnd1
## 188            merNub1
## 189            mesUni1
## 190            micMur1
## 191            micMur2
## 192               mm10
## 193         mm10Patch4
## 194                mm5
## 195                mm6
## 196                mm7
## 197                mm8
## 198                mm9
## 199            monDom1
## 200            monDom4
## 201            monDom5
## 202            musFur1
## 203            myoLuc2
## 204            nanPar1
## 205            nasLar1
## 206            necAme1
## 207            nipNip1
## 208            nomLeu1
## 209            nomLeu2
## 210            nomLeu3
## 211            ochPri2
## 212            ochPri3
## 213            oncVol1
## 214            opiHoa1
## 215            oreNil1
## 216            oreNil2
## 217            oreNil3
## 218            ornAna1
## 219            ornAna2
## 220            oryCun2
## 221            oryLat2
## 222            otoGar3
## 223            oviAri1
## 224            oviAri3
## 225            oviAri4
## 226            panPan1
## 227            panPan2
## 228            panRed1
## 229            panTro1
## 230            panTro2
## 231            panTro3
## 232            panTro4
## 233            panTro5
## 234            panTro6
## 235            papAnu2
## 236            papHam1
## 237            pelCri1
## 238            pelSin1
## 239            petMar1
## 240            petMar2
## 241            petMar3
## 242            phaCar1
## 243            phaLep1
## 244            phoRub1
## 245            picPub1
## 246            ponAbe2
## 247            ponAbe3
## 248            priExs1
## 249            priPac1
## 250            priPac3
## 251            proCap1
## 252     proteins120806
## 253     proteins121210
## 254     proteins140122
## 255     proteins150225
## 256     proteins160229
## 257     proteins180404
## 258           proteome
## 259            pteGut1
## 260            pteVam1
## 261            pygAde1
## 262            pytBiv1
## 263            rheMac1
## 264            rheMac2
## 265            rheMac3
## 266            rheMac8
## 267            rhiRox1
## 268                rn3
## 269                rn4
## 270                rn5
## 271                rn6
## 272            sacCer1
## 273            sacCer2
## 274            sacCer3
## 275            saiBol1
## 276            sarHar1
## 277            serCan1
## 278            sorAra1
## 279            sorAra2
## 280           sp120323
## 281           sp121210
## 282           sp140122
## 283           sp150225
## 284           sp160229
## 285           sp180404
## 286            speTri2
## 287            strCam1
## 288            strPur1
## 289            strPur2
## 290            strRat2
## 291           susScr11
## 292            susScr2
## 293            susScr3
## 294            taeGut1
## 295            taeGut2
## 296            tarSyr1
## 297            tarSyr2
## 298            tauEry1
## 299            tetNig1
## 300            tetNig2
## 301            thaSir1
## 302            tinGut2
## 303            triMan1
## 304            triSpi1
## 305            triSui1
## 306            tupBel1
## 307            turTru2
## 308            tytAlb1
## 309            uniProt
## 310            vicPac1
## 311            vicPac2
## 312           visiGene
## 313            xenLae2
## 314            xenTro1
## 315            xenTro2
## 316            xenTro3
## 317            xenTro7
## 318            xenTro9
## 319            zonAlb1

hg19 <- dbConnect(MySQL(), user="genome", db="hg19",
                  host="genome-mysql.cse.ucsc.edu") #Una base dedatos del genoma humano 19
allTables <- dbListTables(hg19)
length(allTables)# la cantidad total de tablas

## [1] 11142

allTables[1:5]

## [1] "HInv"         "HInvGeneMrna" "acembly"      "acemblyClass"
## [5] "acemblyPep"

dbListFields(hg19,"affyU133Plus2")#las columnas de esta tabla

##  [1] "bin"         "matches"     "misMatches"  "repMatches"  "nCount"     
##  [6] "qNumInsert"  "qBaseInsert" "tNumInsert"  "tBaseInsert" "strand"     
## [11] "qName"       "qSize"       "qStart"      "qEnd"        "tName"      
## [16] "tSize"       "tStart"      "tEnd"        "blockCount"  "blockSizes" 
## [21] "qStarts"     "tStarts"

dbGetQuery(hg19, "select count(*) from affyU133Plus2") # es un comando de Querry para saber cuantas filas tiene la tabla

##   count(*)
## 1    58463

affyData <- dbReadTable(hg19, "affyU133Plus2")# se lee esta tabla y se muetran los primeros datos

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 0 imported
## as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 1 imported
## as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 2 imported
## as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 3 imported
## as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 4 imported
## as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 5 imported
## as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 6 imported
## as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 7 imported
## as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 8 imported
## as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 11
## imported as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 12
## imported as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 13
## imported as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 15
## imported as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 16
## imported as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 17
## imported as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 18
## imported as numeric

head(affyData)

##   bin matches misMatches repMatches nCount qNumInsert qBaseInsert
## 1 585     530          4          0     23          3          41
## 2 585    3355         17          0    109          9          67
## 3 585    4156         14          0     83         16          18
## 4 585    4667          9          0     68         21          42
## 5 585    5180         14          0    167         10          38
## 6 585     468          5          0     14          0           0
##   tNumInsert tBaseInsert strand        qName qSize qStart qEnd tName
## 1          3         898      -  225995_x_at   637      5  603  chr1
## 2          9       11621      -  225035_x_at  3635      0 3548  chr1
## 3          2          93      -  226340_x_at  4318      3 4274  chr1
## 4          3        5743      - 1557034_s_at  4834     48 4834  chr1
## 5          1          29      -    231811_at  5399      0 5399  chr1
## 6          0           0      -    236841_at   487      0  487  chr1
##       tSize tStart  tEnd blockCount
## 1 249250621  14361 15816          5
## 2 249250621  14381 29483         17
## 3 249250621  14399 18745         18
## 4 249250621  14406 24893         23
## 5 249250621  19688 25078         11
## 6 249250621  27542 28029          1
##                                                                   blockSizes
## 1                                                          93,144,229,70,21,
## 2              73,375,71,165,303,360,198,661,201,1,260,250,74,73,98,155,163,
## 3                 690,10,32,33,376,4,5,15,5,11,7,41,277,859,141,51,443,1253,
## 4 99,352,286,24,49,14,6,5,8,149,14,44,98,12,10,355,837,59,8,1500,133,624,58,
## 5                                       131,26,1300,6,4,11,4,7,358,3359,155,
## 6                                                                       487,
##                                                                                                  qStarts
## 1                                                                                    34,132,278,541,611,
## 2                        87,165,540,647,818,1123,1484,1682,2343,2545,2546,2808,3058,3133,3206,3317,3472,
## 3                   44,735,746,779,813,1190,1195,1201,1217,1223,1235,1243,1285,1564,2423,2565,2617,3062,
## 4 0,99,452,739,764,814,829,836,842,851,1001,1016,1061,1160,1173,1184,1540,2381,2441,2450,3951,4103,4728,
## 5                                                     0,132,159,1460,1467,1472,1484,1489,1497,1856,5244,
## 6                                                                                                     0,
##                                                                                                                                      tStarts
## 1                                                                                                             14361,14454,14599,14968,15795,
## 2                                     14381,14454,14969,15075,15240,15543,15903,16104,16853,17054,17232,17492,17914,17988,18267,24736,29320,
## 3                               14399,15089,15099,15131,15164,15540,15544,15549,15564,15569,15580,15587,15628,15906,16857,16998,17049,17492,
## 4 14406,20227,20579,20865,20889,20938,20952,20958,20963,20971,21120,21134,21178,21276,21288,21298,21653,22492,22551,22559,24059,24211,24835,
## 5                                                                         19688,19819,19845,21145,21151,21155,21166,21170,21177,21535,24923,
## 6                                                                                                                                     27542,

query <- dbSendQuery(hg19, "select * from affyU133Plus2 where misMatches between 1 and 3")

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 0 imported
## as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 1 imported
## as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 2 imported
## as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 3 imported
## as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 4 imported
## as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 5 imported
## as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 6 imported
## as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 7 imported
## as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 8 imported
## as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 11
## imported as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 12
## imported as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 13
## imported as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 15
## imported as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 16
## imported as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 17
## imported as numeric

## Warning in .local(conn, statement, ...): Unsigned INTEGER in col 18
## imported as numeric

affyMis <- fetch(query)
quantile(affyMis$misMatches)

##   0%  25%  50%  75% 100% 
##    1    1    2    2    3

affyMisSmall <- fetch(query, n=10)
dbClearResult(query)

## [1] TRUE

dim(affyMisSmall)

## [1] 10 22

dbDisconnect(hg19)

## [1] TRUE

10.2 Leer datos y escribir con HDF5

10.3 Leer datos de páginas web (el código HTML)

10.4 Leer datos de API’s (Application programmming Interface, ejemplo Twitter)

10.5 Leer otro tipo de datos (imagenes, musica, archivos de otrso programs: SAS, Excel)

11 Monte Carlo

Seleccionar muestras aleatorias de una población utilizando números aleatorios

Resolver problemas matemáticos usando muestreo estadístico

Finales de los cuarentas principios de los cincuentas

Stanislau Ulam, En 1946 da el nombre de Monte Carlo en honor de un familiar que es propenso a apostar Cita. wikipedia, enciclopedia Britanica, MacTutor History of Mathematics archive, Mathematics Genealogy Project,
Enrico Fermi, wikipedia, MacTutor History of Mathematics archive, Mathematics Genealogy Project
John Von Neumann, wikipedia, enciclopedia Britanica, MacTutor History of Mathematics archive, Mathematics Genealogy Project
Nicolas Metropolis, wikipedia, wolfram, The New York Times

Para leer

THE BEGINNING of the MONTE CARLO METHOD. Divulgación escrita por Meropolis en 1987
STAN ULAM, JOHN VON NEUMANN, and the MONTE CARLO METHOD. Dibulgación escrita por Roger Eckhardt en 1987
Monte Carlo (Wikipedia)
THE MONTE CARLO METHOD. Investigación escrita por Metropolis y Ulam 1949

11.1 Probabilidad

Probabilidad 1

Probabilidad 3

Probabilidad 4

Probabilidad 5

Probabilidad 6

Probabilidad 7

Probabilidad 8

Probabilidad 9

Probabilidad 10

Probabilidad 11

Probabilidad 12

curve(dunif(x), from=-1, to=2, col="violet",lwd=2,
      xlab="x",ylab="f(x)",main="Función de densidad uniforme")

curve(punif(x), from=-1, to=2, col="orange1", lwd=2,
      xlab="x",ylab="f(x)",main="Función de distribución uniforme")

datos=runif(1000)

plot(datos, xlab="numéro de dato", ylab = "valor del dato",
     col="blue")

boxplot(datos,col="gold",ylab="Rango de valores")

hist(datos, prob = F, col="lightsalmon",main="Histograma",
     sub="Datos simulados de una unif(0,1)",ylab="f(x)",
     breaks = "FD")
curve(dunif(x), from=-1, to=2, add=T,col="violet",lwd=6)

curve(dnorm(x), from=-3, to=3,  col="violet",lwd=2,
      xlab="x",ylab="f(x)",main="Función de densidad normal")

curve(pnorm(x), from=-3, to=3, col="orange1", lwd=2,
      xlab="x",ylab="f(x)",main="Función de distribución normal")

datos=rnorm(1000)

plot(datos, xlab="numéro de dato", ylab = "valor del dato",
     col="blue",pch = 8,cex = 0.1)

boxplot(datos,col="gold",ylab="Rango de valores")

hist(datos, prob = T, col="lightsalmon",main="Histograma",
     sub="Datos simulados de una N(0,1)",ylab="f(x)",
     breaks = "FD")

curve(dnorm(x), from=-3, to=3, add=T,col="violet",lwd=6,
      ylab="f(x)")

curve(dexp(x), from=-1, to=3,  col="violet",lwd=2,
      xlab="x",ylab="f(x)",main="Función de densidad exponencial")

curve(pexp(x), from=-1, to=3, col="orange1", lwd=2,
      xlab="x",ylab="f(x)",main="Función de distribución exponencial")

datos=rexp(1000)

plot(datos, xlab="numéro de dato", ylab = "valor del dato",
     col="blue",pch = "." )

boxplot(datos,col="gold",ylab="Rango de valores")

hist(datos, prob = T, col="lightsalmon",main="Histograma",
     sub="Datos simulados de una exp(1)",ylab="f(x)",
      breaks = "FD")

curve(dexp(x), from=-1, to=5, add=T,col="violet",lwd=6,
      ylab="f(x)")

# Ayuda para distribuciones de probabilidad
# help(Distributions)

Probabilidad 13 Probabilidad 14

Ley de los grandes numéros

leynumeros = function(n)
{
  x=runif(n)
  rango=seq(from=1, to=n, by=1)
  s=c()
  for (i in rango) {
    s[i]=sum(x[1:i])/i
    }
  plot(s, type="l",xlab="n",ylab="S/n",main="Ley de los grandes numéros",col="lightsalmon",sub = "Datos simulados de una distribución unif(0,1)")
}
leynumeros(100)

leynumeros(1000)

leynumeros(10000)

leynumeros(100000)

11.2 Estadística

$Estadistica 1$

Estadistica 1

11.3 Procesos Estocásticos

Procesos Estocásticos 1 Procesos Estocásticos 2 Procesos Estocásticos 3 Procesos Estocásticos 4 Procesos Estocásticos 5

11.4 Monte Carlo

Monte Carlo 1

montecarlo1=function(g,a,b,n)
{
  x=runif(n)
  rango=seq(from=1, to=n, by=1)
  #valoresg=apply(,g)
  valoresg=g((b-a)*x+a)
  s=c()
  for (i in rango) {
    s[i]=(b-a)*sum(valoresg[1:i])/i
  }
  s
  plot(s, type="l",xlab="n",ylab="Estimación integral",main="Monte Carlo para integrales de una variable",col="lightsalmon")
  
numerico=integrate(g, lower = a, upper = b)
lines(c(1,n),c(numerico[1],numerico[1]),col="green",lwd=2)
}

g=function(x)
{
  return((1/sqrt(2*pi))*exp(-(x^2)/2))
}

#montecarlo1(g,-2.5,2.5,10000)

11.5 Generar Variables aleatorias

Monte Carlo 2

Ejemplo: Generar una variable aleatoria con distribución acumulada $F(x)=1-e^{-2 x}$

Solución: \[F^{-1}(x)=-\frac{1}{2}ln(1-x)\]

montecarlo2=function(Fi,n)
{
u=runif(n)
X=Fi(u)
plot(X, xlab="numéro de dato", ylab = "valor del dato",
     col="blue",pch = "." )

hist(X, prob = T, col="lightsalmon",main="Histograma",
     sub="Datos simulados",ylab="f(x)",
      breaks = "FD")

curve(dexp(x,rate=2), from=-1, to=5, add=T,col="violet",lwd=6,
      ylab="f(x)")
}

g=function(x)
{
  return(-(1/2)*log(1-x))
}

montecarlo2(g,10000)

Ejemplo: Generar una variable aleatoria con función de densidad

\[f(x)=\frac{2}{\pi R^2}\sqrt{R^2-x^2},\;\;-R\leq x\leq R\],

con $R=2$

Solución: Graficamos la función de densidad

g=function(x,R)
{
  return((2/(pi*R^2)*(sqrt(R^2-x^2))))
}

grafica=function(R,g)
{
x=seq(from=-R, to=R, by=R/100)
plot(x,g(x,R),col="blue", type = "l", ylab = "f(x)")
}
grafica(2,g)

Se genera $U_1$ y $U_2$ uniformes es el intevalo $(0,1)$ Si $U_1\leq \sqrt{1-(2U_2-1)^2}$ entonces $Z=X=4U_2-2$ si no repetimos.

montecarlo3=function(Fi,n,g)
{
Z=c()
while (length(Z)<n)
{
u1=runif(1)
u2=runif(1);
if (u1<=Fi(u2)){
Z=c(Z,4*u2-2)  
}
}

plot(Z, xlab="numéro de dato", ylab = "valor del dato",
     col="blue",pch = "." )

hist(Z, prob = T, col="lightsalmon",main="Histograma",
     sub="Datos simulados",ylab="f(z)",
      breaks = "FD")
curve(g(x), from=-2, to=2, add=T,col="violet",lwd=6,
      ylab="f(x)")
}

funci=function(x)
{
return(sqrt(1-(2*x-1)^2))
}

g=function(x)
{
  return((2/(pi*2^2)*(sqrt(2^2-x^2))))
}
montecarlo3(funci,10000,g)

11.6 Generar Cadenas de Markov con método de Monte Carlo (MCMC)

El primer algoritmo fue propuesto pot Metropolis (1953), posteriormente Hastings (1970) lo generaliza así el algoritmo se conoce como Metropolis-Hostings (MH) y es el del tipo (MCMC). El algoritmo MH simula la cadena de Markov (calores de $X_t$) cuya distribución estacionaria es $\pi$ (la f chiquita de $X_t$) lo que significa que para $t\rightarrow\infty$ los calores de $X_t$ serán como las muestars de $\pi$.

target = function(x){
  return(ifelse(x<0,0,exp(-x)))
}

x = rep(0,10000)
x[1] = 3     #inicializamos 3
for(i in 2:10000){
  current_x = x[i-1]
  proposed_x = current_x + rnorm(1,mean=0,sd=1)
  A = target(proposed_x)/target(current_x) 
  if(runif(1)<A){
    x[i] = proposed_x       
  } else {
    x[i] = current_x        
  }
}

plot(x,main="valores de x generados con el algoritmo")

hist(x,xlim=c(0,10),probability = TRUE, main="Histograma de los valores de X")
xx = seq(0,10,length=100)
lines(xx,target(xx),col="red")

12 Bibliografía

12.1 Cursos MOOC con R

Statistics and R

Introduction to R for Data Science

Data Science: R Basics

Explorer Statistics with R

Foundations of Data Analysis - Part 1: Statistics Using R

Programming with R for Data Science

Analyzing Big Data with Microsoft R

Statistical Learning

Staistical Modelind and Regression Analysis

12.2 Bibliografia R

Part 1: Beginner’s guide to R: Introduction

Part 2: Getting your data into R

R Markdown Cheat Sheet

R Markdown Reference Guide

Introducción a R

A (very) short introduction to R

R Tutorial

Data Visualization ggplot2

https://rpubs.com/williamsurles

Estadística con R

Quick-R

MATLAB-PYTHON-R

Creating Dynamic Documents with RMarkdown and Knitr

Aprendizaje del Software Estadístico R: un entorno para simulación y computación estadística

El arte de programar en R

12.3 Bibliografía probabilidad

Aunque existe bastante bibliografía al respecto, considero que las notas de probabilidad de mi profesor el Dr. Luis Rincón de la Facultad de ciencias de la UNAM son adecuadas, en el siguiente enlace se puede consultar el libro y este otro enlace los videos que desarrollo el Dr. Luis Rincón sobre el tema.

12.4 Bibliografía estadística

Aunque existe bastante bibliografía al respecto, considero que las notas de estadística de mi profesor el Dr. Luis Rincón de la Facultad de ciencias de la UNAM son adecuadas, en el siguiente se pueden enlace los videos que desarrollo el Dr. Luis Rincón sobre el tema.

12.5 Bibliografía procesos estocásticos

documentos con knit

12.6 Bibliografia Monte Carlo

12.7 Código de las notas

Enlace, las imagenes pueden descargarse de esta presentación para ser utilizadas.

R introducción para Métodos de Monte Carlo

Rafael Martínez Martínez

Marzo, 2018