Ejercicio 1
Utilizando la información correspondiente a los sismos registrados por el Sistema Sismológico Nacional (SSN) en el periodo 01-01-1980 al 24-09-2021, disponible en (http://www2.ssn.unam.mx:8080/catalogo/.html), calcular:
Comenzamos llamando a la base correspondiente:
arreglado<-read.csv("/Users/mac/Downloads/SSNMX_catalogo_19800101_20210924 arreglado.csv")A)
El total de sísmos detectados por mes.
Nota: Por alguna razón el documento .csv con la base de sismos contenía fechas con el formato " 01/01/80" de modo que funciones como as.date no resultaron útiles para la extracción de información. Por lo tanto, (después de algunas horas de frustración y desesperanza) se optó por la función dmy de la paquetería “lubridate” para obtener en una nueva columna, un formato más amigable de manipular.
library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionarreglado$fecha5<- dmy(arreglado$Fecha)
head(arreglado) ## Fecha Hora Magnitud Latitud Longitud Profundidad
## 1 01/01/1980 04:12:00 no calculable 16.98 -102.33 33
## 2 01/01/1980 06:02:37 no calculable 13.50 -91.36 33
## 3 01/01/1980 07:47:55 4.6 13.21 -90.46 59
## 4 01/01/1980 12:30:10 no calculable 16.37 -97.85 40
## 5 01/01/1980 21:42:37 no calculable 16.85 -98.52 33
## 6 02/01/1980 11:16:17 no calculable 16.14 -97.94 77
## Referencia.de.localizacion Fecha.UTC Hora.UTC Estatus
## 1 109 km al SUR de CD LAZARO CARDENAS, MICH 01/01/1980 10:12:00 revisado
## 2 156 km al SURESTE de CD HIDALGO, CHIS 01/01/1980 12:02:37 revisado
## 3 244 km al SURESTE de CD HIDALGO, CHIS 01/01/1980 13:47:55 revisado
## 4 22 km al ESTE de PINOTEPA NACIONAL, OAX 01/01/1980 18:30:10 revisado
## 5 22 km al NOROESTE de OMETEPEC, GRO 02/01/1980 03:42:37 revisado
## 6 25 km al SURESTE de PINOTEPA NACIONAL, OAX 02/01/1980 17:16:17 revisado
## fecha5
## 1 1980-01-01
## 2 1980-01-01
## 3 1980-01-01
## 4 1980-01-01
## 5 1980-01-01
## 6 1980-01-02tail(arreglado)## Fecha Hora Magnitud Latitud Longitud Profundidad
## 226759 24/09/2021 22:27:00 3.5 17.98 -100.03 53
## 226760 24/09/2021 22:59:56 3.7 30.80 -114.34 2
## 226761 24/09/2021 23:05:36 3.7 18.17 -103.15 4
## 226762 24/09/2021 23:09:38 3.5 19.40 -103.93 29
## 226763 24/09/2021 23:28:29 3.5 31.43 -115.52 10
## 226764 24/09/2021 23:47:22 3.7 14.64 -92.20 82
## Referencia.de.localizacion Fecha.UTC Hora.UTC
## 226759 45 km al SURESTE de ARCELIA, GRO 25/09/2021 03:27:00
## 226760 54 km al SURESTE de SAN FELIPE, BC 25/09/2021 03:59:56
## 226761 67 km al SUR de COALCOMAN, MICH 25/09/2021 04:05:36
## 226762 25 km al NOROESTE de VILLA DE ALVAREZ, COL 25/09/2021 04:09:38
## 226763 79 km al NOROESTE de SAN FELIPE, BC 25/09/2021 04:28:29
## 226764 6 km al SUROESTE de CD HIDALGO, CHIS 25/09/2021 04:47:22
## Estatus fecha5
## 226759 verificado 2021-09-24
## 226760 verificado 2021-09-24
## 226761 verificado 2021-09-24
## 226762 verificado 2021-09-24
## 226763 verificado 2021-09-24
## 226764 verificado 2021-09-24Como necesitamos el acumulado de casos por mes, entonces crearemos una nueva columna y usaremos “month” en la columna “fecha5” que recientemente creamos. A continuación, para representar el total de casos por mes, usaremos la función “summarise” de la paquetería “dplyr” que nos resumirá y agrupará los casos en la columna “Mes” que previamente creamos.
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionarreglado$Mes<-lubridate::month(arreglado$fecha5, label=TRUE)
sismos_mens<- summarise(group_by(arreglado, Mes),"Total de sismos" = n())
pander::pander( sismos_mens, style='rmarkdown' ) # style is pipe tables...| Mes | Total de sismos |
|---|---|
| Jan | 21195 |
| Feb | 20644 |
| Mar | 19155 |
| Apr | 17959 |
| May | 17807 |
| Jun | 18129 |
| Jul | 17633 |
| Aug | 17827 |
| Sep | 21497 |
| Oct | 18254 |
| Nov | 18094 |
| Dec | 18570 |
B)
El total de sismos detectados por cada estado de la república.
En este caso necesitamos agrupar a los casos por estado, así que similar al paso anterior, nos será útil crear una columna que identifique la procedencia de cada caso. Para ello, extraeremos los útlimos caracteres de la columna " Referencia de localización".
library(tm)## Loading required package: NLPcolnames(arreglado)## [1] "Fecha" "Hora"
## [3] "Magnitud" "Latitud"
## [5] "Longitud" "Profundidad"
## [7] "Referencia.de.localizacion" "Fecha.UTC"
## [9] "Hora.UTC" "Estatus"
## [11] "fecha5" "Mes"arreglado$estado<-sub(".*,","",arreglado$Referencia.de.localizacion )
sismos_estados <- arreglado %>% group_by(estado) %>%
summarise("Total de sismos" = n())
pander::pander( sismos_estados, style='rmarkdown', caption="Estados" ) | estado | Total de sismos |
|---|---|
| AGS | 42 |
| BC | 7517 |
| BC | 653 |
| BCS | 2784 |
| BCS | 62 |
| CAMP | 62 |
| CDMX | 280 |
| CDMX | 5 |
| CHIH | 399 |
| CHIH | 1 |
| CHIS | 42095 |
| CHIS | 818 |
| COAH | 95 |
| COAH | 14 |
| COL | 5037 |
| COL | 251 |
| DGO | 60 |
| DGO | 1 |
| GRO | 38322 |
| GRO | 927 |
| GTO | 54 |
| GTO | 28 |
| HGO | 536 |
| HGO | 12 |
| JAL | 6427 |
| JAL | 405 |
| MEX | 573 |
| MEX | 10 |
| MICH | 13097 |
| MICH | 882 |
| MOR | 229 |
| MOR | 4 |
| N | 3 |
| NAY | 227 |
| NAY | 36 |
| NL | 405 |
| NL | 15 |
| OAX | 92698 |
| OAX | 1894 |
| PUE | 1105 |
| PUE | 25 |
| QR | 86 |
| QRO | 14 |
| QRO | 1 |
| SIN | 903 |
| SIN | 17 |
| SLP | 178 |
| SLP | 11 |
| SON | 1182 |
| SON | 127 |
| TAB | 508 |
| TAB | 22 |
| TAMS | 112 |
| TAMS | 3 |
| TLAX | 118 |
| TLAX | 2 |
| VER | 5036 |
| VER | 123 |
| YUC | 3 |
| ZAC | 205 |
| ZAC | 23 |
C)
Gráfico circular de frecuencia de los sísmos de acuerdo a la hora del día (0-23 o 1-24) cuando acontenció.
time <-strptime (arreglado$Hora, format = "%H:%M:%S")
arreglado$hr <- hour(time)
arreglado$unos <-1
var_unos<-arreglado$unos
hora<- (aggregate(unos~hr, data = arreglado, sum ))
#hora_n<- as.numeric(hora)library(ggplot2)##
## Attaching package: 'ggplot2'## The following object is masked from 'package:NLP':
##
## annotateggplot(hora, aes(x=hr, y=unos)) +
geom_bar(stat = "identity", fill="#3498DB", position=position_dodge()) +
coord_polar(start = 0)+
scale_x_continuous(breaks = c(0:23), labels = as.character(c(0:23))) +
labs(y= "Casos",
x= "Hora")
D)
Considerando solamente aquellos sismos cuya magnitud fue mayor o igual a 5 ¿cuáles son los 5 estados con mayor cantidad de sısmos?
arreglado$cinco_mas <-as.numeric(arreglado$Magnitud>=5 )
cinco_mas_<-(arreglado$Magnitud>=5)
datos<-(aggregate(cinco_mas ~ estado, arreglado, sum))
final_datos<-datos %>%
top_n(n=5,cinco_mas)%>%
arrange(desc(cinco_mas))
pander::pander( final_datos, style='rmarkdown', caption="Estados con sismos de mayor magnitud") | estado | cinco_mas |
|---|---|
| CHIS | 5695 |
| OAX | 3847 |
| GRO | 2174 |
| VER | 906 |
| MICH | 737 |
E)
- Realizar un histograma de la magnitud de los sismos detectados en el estado de Chiapas.
library(ggplot2)
Magnitud_num<-as.numeric(arreglado$Magnitud)## Warning: NAs introduced by coercionchis_1<- subset(arreglado, estado == "CHIS")
histo<- ggplot(chis_1, aes(x=Magnitud_num, y =var_unos)) +
geom_bar(stat = "identity", fill="#E69F00") + labs(title = "sismos detectados en Chiapas", x= "magnitud", y= "casos")
histo## Warning: Removed 13826 rows containing missing values (position_stack).
Ejercicio 2
El método más común para la determinación de la aproximación estadísticamente òptima con un conjunto de parámetros dado un grupo de datos, es el método conocido como Mínimos Cuadrados (por sus siglas en inglés, LS), y fue propuesto hace aproximadamente dos siglos por Carl Friedrich Gauss. Así, en el caso más simple de la regresión lineal \[(y = b0 + b1X )\] , se busca minimizar la suma de cuadrados en el vector de residuales e. Esta función objetivo de mínimos cuadrados puede ser escrita en su forma compacta de la siguiente manera:
\(f (b) = e^T e\)
\(f (b) = (y - Xb)^T (y-Xb)\)
\(f (b) = y^T y-2y^T Xb + bX^T Xb\)
Calculando las derivadas parciales respecto a los parámetros \(b\), e igualando dicha derivada a 0, se obtiene el vector de estimaciones puntuales de los parámetros \(β\) (ó b):
\[ b= (X^T X)^-1 X^Ty\]
Utilizando los datos disponibles en Moodle llamados , obtener el vector de estimaciones puntuales para los parámetros b cuando X = matriz de variables biking y smoking (498 x 2) y, y = heart.disease (498 x 1).
Comenzamod llamando la base que necesitamos.
corazon<-read.csv("/Users/mac/Downloads/corazon.csv")
colnames(corazon)## [1] "biking" "smoking" "heart.disease"Para encontrar el vector de estimaciones puntuales de los parámetros\[b\] manipularemos a la matriz “corazon” de la siguiente manera:
\[ b= (X^T X)^-1 X^Ty\] La primera parte de la expresión nos pide obtener la inversa de $ (X^T X)$ :
library(matlib)
corazon_m <- as.matrix(corazon)
## creamos una matriz solamente con y (heart disease)
y<- as.matrix(corazon_m[,3])
## modificamos nuestra matriz con variables predictoras para poder tratarla después.
corazon_x<- corazon_m[,-3]Una vez listas, podemos comenzar a tratar nuestras matrices:
library(matlib)
trans<- t(corazon_x)
inv_m<-inv( trans%*% corazon_x )En la segunda parte la expresión se requiere de \(x\) transpuesta que acabamos de obtener y el arreglo con la variable predictora \(y\). Por lo tanto, en conjunto tendríamos:
library(pander)
coef<-inv_m%*% trans%*%y
pander::pander( coef, style='rmarkdown', caption="Parámetros b" ) # style is pipe tables...| -0.03791 |
| 0.6234 |