Integration scanstat.Rmd code

Comparaison_of_methods.rmd

@@ -11,6 +11,7 @@ library("localScore")
 library("latex2exp")
 library("Rcpp")
 library("caret")
+library("ROCR")
 ```
 
 ## 1. Proposition for simulations under $\mathcal{H}_1$
@@ -370,7 +371,7 @@ LocalScoreMC <- function(lambda0, lambda1, NbSeq, T, X_seq, P_X, tbe0){
       
       pvalue = c(pvalue, daudin_result)
   }
-  LS_H0=data.frame(num=1:NbSeq, pvalue_scan=pvalue, class=(pvalue<0.05))
+  LS_H0=data.frame(num=1:NbSeq, pvalue_scan=pvalue, class=as.numeric(pvalue<0.05))
   return(LS_H0)
 }
 ```
@@ -554,3 +555,141 @@ for (i in (1:NbSeq)) {
   titleSpec=TeX(paste(r'(Specificity for $\lambda_0=$)', lambda0))
   plot(x=accepted_lambda,y=Specificity, type='l', main = titleSpec)
 ```
+
+
+
+
+
+## CompareMethods
+```{r}
+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)
+        pp0 =  vector("list", length = NbSeq)
+        pp1 =  vector("list", length = NbSeq)
+        tbe1 = vector("list", length =  NbSeq)
+        
+        for (i in (1:NbSeq)) {
+            #Simulation for sequences under H0
+            ppi = PoissonProcess(lambda0,T)
+            ni=length(ppi)
+            pp0[[i]] = ppi
+            tbei = ppi[2:ni]-ppi[1:ni-1]
+            tbe0[[i]] = tbei
+            
+            #Simulation for sequences under H1
+            ppj1 = SimulationH1(lambda0, lambda1, T, tau)
+            nj = length(ppj1)
+            pp1[[i]] = ppj1
+            tbej = ppj1[2:nj]-ppj1[1:nj-1]
+            tbe1[[i]] = tbej
+        }
+        
+        #cat("- Empiric version:\n")
+        Score = ScoreDistribEmpiric(lambda0, lambda1, NbSeq, 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)
+        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)
+          
+                    
+        cat("--- Confusion matrix for scan statistic method --- \n")
+        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)
+          
+        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)
+        cat("-----------------------------------\n")
+        
+        
+        result <- c('performance.SS'= perf.SS,'performance.LS'= perf.LS)
+        return(result)
+    }
+}
+```
+
+```{r}
+NbSeq = 100
+T = 10
+tau = 2
+lambda0 = 0.3
+lambda1 = 0.5
+
+result1 = CompareMethods(lambda0, lambda1, NbSeq, T, tau)
+
+lambda0 = 0.01
+lambda1 = 1
+
+result2 = CompareMethods(lambda0, lambda1, NbSeq, T, tau)
+
+lambda0 = 1
+lambda1 = 1.1
+
+result3 = CompareMethods(lambda0, lambda1, NbSeq, T, tau)
+
+
+lambda0 = 0.9
+lambda1 = 2
+
+result4 = CompareMethods(lambda0, lambda1, NbSeq, T, tau)
+
+
+title_ROC = TeX(paste(r'(ROC curve for several values of $\lambda_0$ and $\lambda_1$)'))
+
+perf1SS = result1[1]
+perf1LS = result1[2]
+                
+perf2SS = result2[1]
+perf2LS = result2[2]
+
+perf3SS = result3[1]
+perf3LS = result3[2]
+
+perf4SS = result4[1]
+perf4LS = result4[2]
+
+
+plot(perf1SS$performance.SS, lty=1, col="coral", lwd = 2)
+par(new=T)
+plot(perf1LS$performance.LS, lty=2, col="coral",lwd = 2)
+
+par(new=T)
+plot(perf2SS$performance.SS, lty=1, col="cyan4", lwd = 2)
+par(new=T)
+plot(perf2LS$performance.LS, lty=2, col="cyan4", lwd = 2)
+
+par(new=T)
+plot(perf3SS$performance.SS, lty=1, col="magenta4", lwd = 2)
+par(new=T)
+plot(perf3LS$performance.LS, lty=2, col="magenta4", lwd = 2)
+
+par(new=T)
+plot(perf4SS$performance.SS, lty=1, col="olivedrab4", lwd = 2)
+par(new=T)
+plot(perf4LS$performance.LS, lty=2, col="olivedrab4", lwd = 2,main=title_ROC)
+
+legend(0.5, 0.3, legend=c("lambda0 = 0.3, lambda1 = 0.5", "lambda0 = 1, lambda1 = 9", "lambda0 = 2, lambda1 = 6", "lambda0 = 8, lambda1 = 9", "lambda0 = 0.1, lambda1 = 0.2")
+       ,col=c("coral", "cyan4", "magenta4", "olivedrab4", "lightgoldenrod3"),  lty=1, cex=0.9,lwd=4,
+       box.lty=0)
+```