From 864fd506a1fd0107eb5ae283b95e6f54df2e9d4a Mon Sep 17 00:00:00 2001
From: Paul-Corbalan <58653590+Paul-Corbalan@users.noreply.github.com>
Date: Tue, 10 May 2022 11:21:14 +0200
Subject: [PATCH] Integration of tau in compare and PValue debug

---
 Comparaison_of_methods.rmd | 39 +++++++++++++++++++++-----------------
 1 file changed, 22 insertions(+), 17 deletions(-)

diff --git a/Comparaison_of_methods.rmd b/Comparaison_of_methods.rmd
index 8ccab51..92e04d0 100644
--- a/Comparaison_of_methods.rmd
+++ b/Comparaison_of_methods.rmd
@@ -12,6 +12,7 @@ library("latex2exp")
 library("Rcpp")
 library("caret")
 library("ROCR")
+library("pROC")
 ```
 
 ## 1. Proposition for simulations under $\mathcal{H}_1$
@@ -120,17 +121,18 @@ CDF=Plot_CDF(lambda0,n_sample,T,tau)
 Compute $p$-value for scan statistic of `ppH1` with `Emp`:
 ```{r}
 PValue <- function(Emp,ppH, T, tau){
-    scanH1=ScanStat(ppH,T,tau)[2]
-    index_scanH1=ScanStat(ppH,T,tau)[1]
-    index=Emp$index_scan
+    SS = ScanStat(ppH,T,tau)
+    scanH = SS[2]
+    index_scanH = SS[1]
+    index = Emp$index_scan
     n=length(index)
-    if (scanH1< min(Emp$index_scan)){
-        return (c(scanH1,1,index_scanH1))
+    if (scanH< min(Emp$index_scan)){
+        return (c(scanH,1,index_scanH))
         } else{
-            if(min(Emp$index_scan)<scanH1 && scanH1<=max(Emp$index_scan)){
-                return(c(scanH1,1-Emp$cdf[scanH1-min(Emp$index_scan)+1],index_scanH1))
-            } else{return (c(scanH1,0,index_scanH1))}}
-}
+            if(min(Emp$index_scan)<scanH && scanH<=max(Emp$index_scan)){
+                return(c(scanH,1-Emp$cdf[scanH-min(Emp$index_scan)],index_scanH))
+            } else{return (c(scanH,0,index_scanH))}}
+    }
 ```
 
 ### 2.2. Simulation under $\mathcal{H}_0$ and computation of p-values
@@ -382,7 +384,7 @@ CompareMethods <- function(lambda0, lambda1, NbSeq, T, tau){
     if (lambda0 < lambda1){
         
         cat("For T = ", T, ", Nb = ", NbSeq, ", lambda0 = ", lambda0, " and lambda1 = ", lambda1, ":\n", sep = "")
-        tbe0 = vector("list", length = NbSeq)
+        tbe0 = vector("list",length=NbSeq)
         pp0 =  vector("list", length = NbSeq)
         pp1 =  vector("list", length = NbSeq)
         tbe1 = vector("list", length =  NbSeq)
@@ -396,7 +398,7 @@ CompareMethods <- function(lambda0, lambda1, NbSeq, T, tau){
             tbe0[[i]] = tbei
             
             #Simulation for sequences under H1
-            ppj1 = SimulationH1(lambda0, lambda1, T, tau)
+            ppj1 = SimulationH1(lambda0, lambda1, T, 3)
             nj = length(ppj1)
             pp1[[i]] = ppj1
             tbej = ppj1[2:nj]-ppj1[1:nj-1]
@@ -404,13 +406,13 @@ CompareMethods <- function(lambda0, lambda1, NbSeq, T, tau){
         }
         
         #cat("- Empiric version:\n")
-        Score = ScoreDistribEmpiric(lambda0, lambda1, NbSeq, T)
+        Score = ScoreDistribEmpiric(lambda0, lambda1, 10**4, T)
         LS_H0 = LocalScoreMC(lambda0, lambda1, NbSeq, T, Score$Score_X, Score$P_X, tbe0)
         LS_H1 = LocalScoreMC(lambda0, lambda1, NbSeq, T, Score$Score_X, Score$P_X, tbe1)
         LS_obtained = c(LS_H0$class, LS_H1$class)
         options(warn = -1) 
           
-        Emp = EmpDistrib(lambda0,n_sample,T,tau)
+        Emp = EmpDistrib(lambda0,10**4,T,tau)
         SS_H0 = ScanStatMC(NbSeq, T, tau, Emp, pp0)
         SS_H1 = ScanStatMC(NbSeq, T, tau, Emp, pp1)
         SS_obtained = c(SS_H0$class, SS_H1$class)
@@ -420,24 +422,27 @@ CompareMethods <- function(lambda0, lambda1, NbSeq, T, tau){
         theoretical_results_SS = c(rep(0,length(SS_H0$num)), rep(1,length(SS_H1$num)))
         print(confusionMatrix(as.factor(SS_obtained), as.factor(theoretical_results_SS),
                               dnn = c("Prediction", "Reference"))$table)
+        roc_SS = roc(theoretical_results_SS, SS_obtained)
+        areaSS = auc(roc_SS)
+        cat("Area under the ROC curve for SS = ", areaSS, "\n")
           
         cat("--- Confusion matrix for local score method --- \n")
         theoretical_results_LS = c(rep(0,length(LS_H0$num)), rep(1,length(LS_H1$num)))
         print(confusionMatrix(as.factor(LS_obtained), as.factor(theoretical_results_LS),
                                 dnn = c("Prediction", "Reference"))$table)
           
-        #cat("--- Coube ROC associé")
         title_ROC = TeX(paste(r'(ROC curve for $H_0: \lambda_0=$)', lambda0, 
                                 r'(against $H_1: \lambda_0=$)', lambda1))
         pred.SS = prediction(theoretical_results_SS,SS_obtained)
         pred.LS = prediction(theoretical_results_LS,LS_obtained)
         perf.SS = performance(pred.SS,"tpr", "fpr")
         perf.LS = performance(pred.LS,"tpr", "fpr")
-        #plot(perf.SS, lty=1, col="coral")
         par(new=T)
-        #plot(perf.LS, lty=2, col="coral", main=title_ROC)
+        roc_LS = roc(theoretical_results_LS, LS_obtained)
+        areaLS = auc(roc_LS)
+        cat("Area under the ROC curve for LS = ", areaLS, "\n")
         cat("-----------------------------------\n")
-        
+        options(warn = -1)        
         
         result <- c('performance.SS'= perf.SS,'performance.LS'= perf.LS)
         return(result)