Dataset Ganado Caprino 7/9/18

library(readr)
library(dplyr)
library(lubridate)
library(ggplot2)
library(plotly)
library(ggmap)

rad2deg <- function(rad) {(rad * 180) / (pi)}
deg2rad <- function(deg) {(deg * pi) / (180)}



caprino_hours<-read.csv("data/caprino_hours.csv.gz",header = F)
caprino_imu_gyros<-read.csv("data/caprino_imu_gyro.csv.gz", header=F)
caprino_imu_acc<-read.csv("data/caprino_imu_acc.csv.gz", header=F)
caprino_gps<-read.csv("data/caprino_gps.csv",header=F)

caprino_imu_acc %>% nrow()
caprino_imu_gyros %>% nrow()
caprino_gps %>% nrow()
caprino_hours %>% nrow()

Datos de los archivos

• caprino_imu_acc.csv, 4 columnas: tiempo UTC en microsegundos, aceleración X en m/s², aceleración Y en m/s², aceleración, Z en m/s². Muestreado a 200 muestras por segundo.

• caprino_imu_gyro.csv, 4 columnas: tiempo UTC en microsegundos, velocidad angular X en rad/s, velocidad angular Y en rad/s, velocidad angular Z en rad/s. Muestreado a 200 muestras por segundo.

• caprino_gps.csv, 4 columnas: tiempo UTC en microsegundos, latitud en radianes, longitud en radianes, altitud en metros. Muestreado a 5 muestras por segundo.

• caprino_hours.csv: tiempo en horas, minutos y segundos. Muestreado a 5 muestras por segundo. La primera hora de este archivo coincide con el primer valor de tiempo UTC del archivo caprino_gps.csv.

Construccion del Dataset

caprino_data<-cbind(caprino_imu_acc,caprino_imu_gyros %>% select(-V1))
caprino_data<-cbind(caprino_data[seq(2,nrow(caprino_data),2),] %>% head(nrow(caprino_gps)-1),  caprino_gps[2:nrow(caprino_gps),] ) 
names(caprino_data)<-c("imu_ts",
                       "accX",
                       "accY",
                       "accZ",
                       "gyroX",
                       "gyroY",
                       "gyroZ",
                       "gps_hour",
                       "gps_ts",
                       "lat",
                       "long",
                       "alt"
                        )
#caprino_data %>% mutate(diff=imu_ts-gps_ts) %>% select(diff)
caprino_data<-caprino_data %>% mutate(activity=ifelse(hms(gps_hour)>hms("8:56:00") & hms(gps_hour)<hms("10:38:00"),"Dco",
                                       
                                        ifelse(hms(gps_hour)>=hms("10:38:00") & hms(gps_hour)<hms("11:38:00"),"C",
                                        ifelse(hms(gps_hour)>=hms("11:38:00") & hms(gps_hour)<hms("12:05:00"),"P",
                                        ifelse(hms(gps_hour)>=hms("12:05:00") & hms(gps_hour)<hms("15:48:00"),"Dca",
                                        ifelse(hms(gps_hour)>=hms("17:30:00") & hms(gps_hour)<hms("20:33:00"),"P",
                                        ifelse(hms(gps_hour)>=hms("20:50:00") & hms(gps_hour)<hms("21:31:00"),"C",
                                        ifelse(hms(gps_hour)>=hms("21:33:00") & hms(gps_hour)<hms("21:50:00"),"Dco",NA)
                                        )))))
                                        )
                        
                        )
caprino_data<-caprino_data %>% mutate(id=cut(as.POSIXct(hms(caprino_data$gps_hour),origin="1960-01-01",tz="GMT"),breaks="1 min")) 

Numero de mediciones por tipo de actividad

NA hace referencia a aquellas mediciones que no fueron etiquetadas

caprino_data %>% group_by(activity) %>% summarise(n=n())

Grafica de la diferencia entre el tiempo de muestreo de la IMU y GPS.

  ggplot(caprino_data%>% mutate(diff=imu_ts - gps_ts))+
  geom_line(aes(y=diff,x=seq(1,nrow(caprino_data))),color='blue')+
  theme_bw()

Grafica del trayecto georeferenciado discriminado por actividad

Grafica del trayecto georeferenciado discriminado por actividad

Las zonas con menor densidad de puntos indican un menor lapso de tiempo en esa zona.

caprino_data   %>%
ggplot(aes(x=long,y=lat))+
  geom_point(aes(color=activity),alpha=0.01,size=0.8)+
  scale_color_brewer(palette="Spectral")+
  theme_bw()

caprino_data_altitude_group<-caprino_data %>% mutate(altitude_group = cut(alt, breaks = c(200,378,525,527,530,650,800,900))) 

ggplot(caprino_data_altitude_group,aes(x=long,y=lat))+
  stat_density2d(aes(fill=altitude_group,alpha=..level..), geom="polygon",  bins=1000) +
  theme_minimal()
  

Distribución de los resultados de los sensores en cada una de las actividades (sin outliers)

• Se observan diferencias en la distribucion de la actividad P en comparacion a las otras actividades.
• Las actividades Dca y Dco presentan similtudes (lo cual es esperable).
• La activdad P presenta algunas similitudes con Dca y Dco, lo que puede dificultar su detection.

caprino_data_melted<-caprino_data %>% select(accX,accY,accZ,gyroX,gyroY,gyroZ,activity) %>% reshape2::melt() 

Using activity as id variables

  ggplot(caprino_data_melted)+
  geom_boxplot(aes(x=variable,color=variable,y=value),outlier.size=0.1,outlier.shape = NA)+
  scale_y_continuous(limits = c(-25,10))+
  theme_bw()+
theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  facet_wrap(~activity) 

NA

Tomamos muestras aleatorias de 1 minuto de duracion para sobre las cuales despues querremos predecir el comportamiento

set.seed(21092019)
minute_slice<-as.vector(t(caprino_data %>% group_by(id) %>% summarise(n=n()) %>% select(id)))
trainidx<-sample(1:length(minute_slice),(length(minute_slice)*70)/100.0)
data_train<-caprino_data %>% filter( id %in% minute_slice[trainidx])
data_test<-caprino_data %>% filter( id %in% (minute_slice[-trainidx]))
data_test%>% group_by(id) %>% summarise(n=n())

Aplicamos Bag-of-features.

presupuesto

Consideramos cada muestra como un vector de 5 dimensiones. La idea es aplicar un algoritmo de clustering sobre estos vectores (accX,accY,accZ,GyroX,gyroY). Una vez encontrado los centroides, para cada muestra de 1 se construye un histograma de los vectores que han caido en cada uno de estos clusters. Este histograma es el que se utiliza para entregar el algoritmo.

La hipotesis de aplicar este metodo es que aquellas muestras que pertenecen a una misma actividad tendran valores similares. Los cuales han sido representado mediantes estos histogramas.

Aplicamos k-means con 12 centroides sobre el conjunto de entrenamiento.

cluster_results<-kmeans(x = data_train%>% select(accX,accY,accZ,gyroX,gyroY), centers = 24, nstart = 5)

did not converge in 10 iterationsQuick-TRANSfer stage steps exceeded maximum (= 8090550)Quick-TRANSfer stage steps exceeded maximum (= 8090550)Quick-TRANSfer stage steps exceeded maximum (= 8090550)Quick-TRANSfer stage steps exceeded maximum (= 8090550)

data_train_cluster<-cbind(data_train,cluster=cluster_results$cluster) %>% group_by(activity,cluster)
data_train_cluster %>% summarise(n=n()) %>%
  ggplot()+
  geom_col(aes(x=cluster,y=n),fill='skyblue')+
  facet_wrap(~activity)+
  theme_bw()

NA

Aplicamos k-means sobre el conjunto de test. Se usan los mismos centroides encontrados en training

cluster_results_test<-kmeans(x = data_test%>% select(accX,accY,accZ,gyroX,gyroY), centers = cluster_results$centers)
data_test_cluster<-cbind(data_test,cluster=cluster_results_test$cluster) %>% group_by(activity,cluster)
data_test_cluster %>% summarise(n=n()) %>%
  ggplot()+
  geom_col(aes(x=cluster,y=n),fill='orange')+
  facet_wrap(~activity)+
  theme_bw()

Preparando el dataset para aplicar un clasificador

Tomamos intervalos de 1m y creamos los atributos usando bag-of-features

sample_minute_slice<-minute_slice[sample(1:length(minute_slice),40)]
data_train_cluster %>% filter (id %in% sample_minute_slice)  %>% group_by(activity,id,cluster) %>% summarise(n=n()) %>%
  ggplot()+
  geom_col(aes(x=cluster,y=n,fill=activity))+
  facet_wrap(~id)+
  theme_bw()

data_train_cluster_labeled<-data_train_cluster   %>% group_by(id,cluster) %>% summarise(n=n()) %>% reshape2::dcast(id~cluster,fill = 0) %>% inner_join(data_train_cluster%>% ungroup() %>% select(id,activity)  %>% unique()  ,by="id") %>% filter(!is.na(activity))

Using n as value column: use value.var to override.

data_test_cluster_labeled<-data_test_cluster   %>% group_by(id,cluster) %>% summarise(n=n()) %>% reshape2::dcast(id~cluster,fill = 0) %>% inner_join(data_test_cluster%>% ungroup() %>% select(id,activity)  %>% unique() ,by="id") 

Using n as value column: use value.var to override.

Probamos Random Forest

Aplicamos 2x5 CV para elegir la mejor version del modelo

library(caret)
library(doMC)
registerDoMC(cores=6)
ctrl_fast <- trainControl(method="repeatedcv", 
                     repeats=2,
                     number=5, 
                     #summaryFunction=twoClassSummary,
                     verboseIter=F,
                     #preProcOptions = list(pcaComp = 10),
                     classProbs=TRUE,
                     allowParallel = T)  
nnGrid <-  expand.grid(decay = c(0), size=c(10,20,50,100,200,300)) 
svmGrid <-  expand.grid(sigma = c(0.001,0.00001), C=c(0.1,0.001,0.5)) 
kerasGrid <- expand.grid(
  size=c(20,50,2000,5000),
  lambda=c(0),
  batch_size=c(1024),
  decay=c(0),
  activation=c('tanh','sigmoid','relu'),
  lr=c(0.001), rho=c(0.9)
)
train_formula=formula("activity~.")
#ctrl_fast$sampling<-"smote"
rfFit <- caret::train(train_formula,
               data = data_train_cluster_labeled %>% select(-id),
               #metric="ROC",
               #preProcess=c("pca"),
               
               #method = "mlpWeightDecay",
               #method = "svmRadial",
               #method ="mlpKerasDecay",
               method = "rf",
               #tuneGrid=kerasGrid,
               tuneLength=11,
               verbose=0,
               #epochs=100,
               trControl = ctrl_fast)

modelos obtenidos

rfFit

Random Forest 

459 samples
 24 predictors
  4 classes: 'C', 'Dca', 'Dco', 'P' 

No pre-processing
Resampling: Cross-Validated (5 fold, repeated 2 times) 
Summary of sample sizes: 366, 368, 368, 368, 366, 368, ... 
Resampling results across tuning parameters:

  mtry  Accuracy   Kappa    
   2    0.7840887  0.7014212
   4    0.7951013  0.7175083
   6    0.7983624  0.7220231
   8    0.7907174  0.7116646
  10    0.7896777  0.7099344
  13    0.7830726  0.7005419
  15    0.7853416  0.7040351
  17    0.7831435  0.7006733
  19    0.7841949  0.7021027
  21    0.7831435  0.7010262
  24    0.7754748  0.6901564

Accuracy was used to select the optimal model using the largest value.
The final value used for the model was mtry = 6.

Evaluacion sobre el conjunto de test

Matriz de Confusion

data_test_prediction<-predict(rfFit,data_test_cluster_labeled %>% select(-id))#,type="prob")
caret::confusionMatrix(data_test_prediction,as.factor(data_test_cluster_labeled$activity))

Confusion Matrix and Statistics

          Reference
Prediction  C Dca Dco  P
       C   25   0   1  1
       Dca  0  65   5  9
       Dco  0   1  27  4
       P    1   7   4 44

Overall Statistics
                                          
               Accuracy : 0.8299          
                 95% CI : (0.7695, 0.8799)
    No Information Rate : 0.3763          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.7605          
 Mcnemar's Test P-Value : NA              

Statistics by Class:

                     Class: C Class: Dca Class: Dco Class: P
Sensitivity            0.9615     0.8904     0.7297   0.7586
Specificity            0.9881     0.8843     0.9682   0.9118
Pos Pred Value         0.9259     0.8228     0.8438   0.7857
Neg Pred Value         0.9940     0.9304     0.9383   0.8986
Prevalence             0.1340     0.3763     0.1907   0.2990
Detection Rate         0.1289     0.3351     0.1392   0.2268
Detection Prevalence   0.1392     0.4072     0.1649   0.2887
Balanced Accuracy      0.9748     0.8874     0.8489   0.8352

CONCLUSIONES