library(plyr)
library(Rcmdr)
library(ggplot2)

#generation echantillon sd = 0.5

x1<- data.frame( rnorm(2000,mean=5, sd=0.5))
names(x1)[c(1)] <- c("Bark")

x1$Groupe<-0
x2<-data.frame( rnorm(2000,mean=7, sd=0.5))
names(x2)[c(1)] <- c("Bark")
x2$Groupe<-1
ROC05<-rbind(x1,x2)

ROC05$Comportement <- factor(ROC05$Groupe, labels=c('Fuite','Attaque'))

p<-ggplot(ROC05,aes(x = Bark, group = Comportement, fill = Comportement))+geom_histogram(aes(x=Bark,y=..density..),alpha=0.5,position="identity")+scale_x_continuous(breaks = seq(-5, 15, 2), lim = c(-5, 15))
p<-p+ggtitle(label="Distribution de l'intensité des aboiements en cas d'attaque et de fuite", subtitle = "sd= 0.5")
 p<-p+scale_fill_manual( values=c("brown","red","green"))
 p
 t.test(Bark~Comportement, alternative='two.sided', conf.level=.95,
 var.equal=FALSE, data=ROC05)

p<-p+geom_vline(xintercept=4.98, color = "brown", size = 1)
 p05<-p+geom_vline(xintercept=7.02, color = "red", size =1)+scale_y_continuous(limits = c(0, 0.80))
 p05

#generation echantillon sd = 3

x1<- data.frame( rnorm(2000,mean=5, sd=3))
 names(x1)[c(1)] <- c("Bark")
x1$Groupe<-0
x2<-data.frame( rnorm(2000,mean=7, sd=3))
names(x2)[c(1)] <- c("Bark")
x2$Groupe<-1

ROC30<-rbind(x1,x2)
 ROC30$Comportement <- factor(ROC30$Groupe, labels=c('Fuite','Attaque'))

p<-ggplot(ROC30,aes(x = Bark, group = Comportement, fill = Comportement))+geom_histogram(aes(x=Bark,y=..density..),alpha=0.5,position="identity")+scale_x_continuous(breaks = seq(-5, 15, 2), lim = c(-5, 15))
 p<-p+ggtitle(label="Distribution de l'intensité des aboiements en cas d'attaque de fuite", subtitle = "sd= 3.0")
 p<-p + scale_y_continuous(limits = c(0, 0.80))+ scale_fill_manual( values=c("brown","red","green"))

t.test(Bark~Comportement, alternative='two.sided', conf.level=.95,
 var.equal=FALSE, data=ROC30)

p<-p+geom_vline(xintercept=4.98, color = "brown", size = 1)
 p30<-p+geom_vline(xintercept=7.02, color = "red", size =1)
 p30
 library(gridExtra)
 grid.arrange(p30, p05, ncol=2)

#analyse ROC
 library(ROCR)

par(mfrow=c(1,1))
 pred1 <- prediction(ROC05$Bark,ROC05$Groupe)
 class(pred1)
 roc.perf1 = performance(pred, measure = "tpr", x.measure = "fpr")
 plot(roc.perf1)
 abline(a=0, b= 1)

pred2 <- prediction(ROC30$Bark,ROC30$Groupe)
 class(pred2)
 roc.perf2 = performance(pred2, measure = "tpr", x.measure = "fpr")
 plot(roc.perf2)

heatcols <- heat.colors(9)
 heatcols <- rainbow(11)

plot(roc.perf1, colorize = TRUE,colorkey=TRUE,colorize.palette=heatcols,lwd=5)
 plot(roc.perf2, add = TRUE, colorize = TRUE, colorkey=TRUE,colorize.palette=heatcols, lwd=5)
 abline(a=0, b= 1)

## precision/recall curve (x-axis: recall, y-axis: precision)
 par(mfrow=c(1,2))

perf1 <- performance(pred1, "rec", "prec")
perf2 <- performance(pred2, "rec", "prec")

plot(perf2, colorize = TRUE,colorkey=TRUE,colorize.palette=heatcols,lwd=4)
plot(perf1, add = TRUE, colorize = TRUE, colorkey=TRUE,colorize.palette=heatcols, lwd=4)

## sensitivity/specificity curve (x-axis: specificity,
 ## y-axis: sensitivity)

perf2 <- performance(pred2, "sens", "spec")
 perf1 <- performance(pred1, "sens", "spec")

plot(perf2, colorize = TRUE,colorkey=TRUE,colorize.palette=heatcols,lwd=4)
 plot(perf1, add = TRUE, colorize = TRUE, colorkey=TRUE,colorize.palette=heatcols, lwd=4))

#from https://www.r-bloggers.com/a-small-introduction-to-the-rocr-package/

opt.cut = function(perf, pred){
 cut.ind = mapply(FUN=function(x, y, p){
 d = (x - 0)^2 + (y-1)^2
 ind = which(d == min(d))
 c(sensitivity = y[[ind]], specificity = 1-x[[ind]],
 cutoff = p[[ind]])
 }, perf@x.values, perf@y.values, pred@cutoffs)
 }

print(opt.cut(roc.perf1, pred1))
 print(opt.cut(roc.perf2, pred2))

cost.perf1 = performance(pred1, "cost", cost.fp = 5, cost.fn = 1)
 pred@cutoffs[[1]][which.min(cost.perf1@y.values[[1]])]

cost.perf2 = performance(pred2, "cost", cost.fp = 5, cost.fn = 1)
 pred@cutoffs[[1]][which.min(cost.perf2@y.values[[1]])]