From 40f1d8cba140468975e2d365f5e2c35014bfafe7 Mon Sep 17 00:00:00 2001
From: elisaduz <elisa.duz1807@gmail.com>
Date: Mon, 9 May 2022 17:15:00 +0200
Subject: [PATCH] Create scanstat.Rmd

---
 scanstat.Rmd | 478 +++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 478 insertions(+)
 create mode 100644 scanstat.Rmd

diff --git a/scanstat.Rmd b/scanstat.Rmd
new file mode 100644
index 0000000..1278839
--- /dev/null
+++ b/scanstat.Rmd
@@ -0,0 +1,478 @@
+---
+title: "scanstat"
+output: pdf_document
+date: '2022-05-09'
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+```
+
+```{r}
+library("localScore")
+library("latex2exp")
+library("Rcpp")
+library("caret")
+library("ROCR")
+```
+
+```{r}
+PoissonProcess <- function(lambda,T) {
+  return(sort(runif(rpois(1,lambda*T),0,T)))
+}
+```
+
+```{r}
+SimulationH1 <- function(lambda0, lambda1,T,tau){
+    ppH0=PoissonProcess(lambda0,T)
+    ppH1.segt=PoissonProcess(lambda1,tau)
+    dbt=runif(1,0,T-tau)
+    ppH0bis=PoissonProcess(lambda0,T)
+    ppH1.repo=dbt+ppH1.segt
+    ppH0_avant=ppH0bis[which(ppH0bis<ppH1.repo[1])]
+    ppH0_apres=ppH0bis[which(ppH0bis>ppH1.repo[length(ppH1.repo)])]
+    ppH1=c(ppH0_avant,ppH1.repo,ppH0_apres)
+    return (ppH1)
+}
+```
+
+```{r}
+TimeBetweenEvent <- function(pp){
+    n=length(pp)
+    tbe=pp[2:n]-pp[1:n1-1]
+    tbe=c(0,tbe)
+    return (tbe)
+}
+
+DataFrame <- function(pp,tbe){
+    list=data.frame(ProcessusPoisson=pp, TimeBetweenEvent=tbe)
+}
+```
+
+```{r}
+ScanStat <- function(pp, T, tau){
+    n=length(pp)
+    stop=n-length(which(pp>(T-tau)))
+    ScanStat=0
+    for (i in (1:stop)) {
+        x=which((pp>=pp[i])&(pp<=(pp[i]+tau)))
+        scan=length(x)
+        if (scan>ScanStat) {ScanStat=scan
+        max=i}
+  }   
+    return (c(max,ScanStat))
+}
+```
+
+```{r}
+EmpDistrib <- function(lambda, n_sample,T,tau){
+    pp=PoissonProcess(lambda,T)
+    scan=c(ScanStat(pp,T, tau)[2])
+    index=c(ScanStat(pp,T, tau)[1])
+    for (i in 2:(n_sample)){
+        pp=PoissonProcess(lambda,T)
+        scan=rbind(scan,ScanStat(pp,T, tau)[2])
+        index=rbind(index,ScanStat(pp,T, tau)[1])
+    }
+    min_scan=min(scan)-1
+    max_scan=max(scan)
+    table1=table(factor(scan, levels = min_scan:max_scan))
+    EmpDis=data.frame(cdf=cumsum(table1)/sum(table1), proba=table1/sum(table1), index_scan=min_scan:max_scan)
+    EmpDis<-EmpDis[,-2]
+    return(EmpDis)
+    }
+```
+
+
+```{r}
+Plot_CDF <- function(lambda,n_sample,T,tau){
+    Emp=EmpDistrib(lambda,n_sample,T,tau)
+    title=TeX(paste(r'(Cumulative distribution function for $\lambda=$)', lambda))
+    plot(Emp$index_scan, Emp$cdf,type="s",xlab="Number of occurrences",ylab="Probability", main=title, col="red")
+    return(Emp)
+}
+```
+
+```{r}
+n_sample=10**4
+lambda0=3
+T=10
+tau=1
+ppH0=PoissonProcess(lambda0,T)
+#CDF=Plot_CDF(lambda0,n_sample,T,tau)
+```
+
+```{r}
+PValue <- function(Emp,ppH, T, tau){
+    scanH1=ScanStat(ppH,T,tau)[2]
+    index_scanH1=ScanStat(ppH,T,tau)[1]
+    index=Emp$index_scan
+    n=length(index)
+    if (scanH1< min(Emp$index_scan)){
+        return (c(scanH1,1,index_scanH1))
+        } 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))}}
+    }
+```
+
+```{r}
+NbSeqH0 = 10
+NbSeqH1 = NbSeqH0
+DataH0 = vector("list")
+DataH1 = vector("list")
+lambda0 = 2
+lambda1 = 5
+T = 10
+tau = 1
+
+#Creation of a sequence that contains the sequence simulated under the null hypothesis
+for (i in 1:NbSeqH0) {
+    ppi = PoissonProcess(lambda0,T)
+    DataH0[[i]] = ppi
+}
+
+#Creation of a sequence that contains the sequence simulated under the alternative hypothesis
+seqH1begin = c()
+for (i in 1:NbSeqH1) {
+    pphi = SimulationH1(lambda0, lambda1,T,tau)
+    DataH1[[i]] = pphi
+}
+
+#Computation of the time between events
+TimeBetweenEventList <- function(list,n_list){
+    TBE = vector("list",length=n_list)
+    for (i in (1:n_list)) {
+        ppi = list[[i]]
+        ni = length(ppi)
+        tbei = ppi[2:ni]-ppi[1:ni-1]
+        TBE[[i]] = tbei
+    }
+    return (TBE)
+}
+tbe0 = TimeBetweenEventList(DataH0,NbSeqH0)
+```
+
+```{r}
+#We start by computing the empirical distribution for lambda0
+Emp = EmpDistrib(lambda0,n_sample,T,tau)
+scan = c()
+pvalue = c()
+index_scan = c()
+
+#Then, we stock the p-value and the 
+for (i in 1:NbSeqH0){
+    ppi = DataH0[[i]]
+    result = PValue(Emp,ppi,T,tau)
+    scan = c(scan,result[1])
+    pvalue = c(pvalue,result[2])
+    index_scan = c(index_scan,result[3])
+}
+
+ScS_H0=data.frame(num=(1:NbSeqH0), scan_stat=scan, pvalue_scan=pvalue,class=c(pvalue<0.05)*1) 
+head(ScS_H0)
+sum(ScS_H0$class[which(ScS_H0$class=='1')])/NbSeqH0
+```
+
+```{r}
+#We start by computing the empirical distribution for lambda0
+scan=c()
+pvalue=c()
+index_scan=c()
+
+#Then, we stock the p-value and the 
+for (i in 1:NbSeqH1){
+    ppi=DataH1[[i]]
+    result=PValue(Emp,DataH1[[i]],T,tau)
+    scan=c(scan,result[1])
+    pvalue=c(pvalue,result[2])
+    index_scan=c(index_scan,result[3])
+}
+ScS_H1 = data.frame(num=1:NbSeqH1, scan_stat=scan, pvalue_scan=pvalue, class=(pvalue<0.05)*1, begin_scan=index_scan)
+head(ScS_H1)
+sum(ScS_H1$class[which(ScS_H1$class=='1')])/NbSeqH1
+```
+
+```{r}
+ScanStatMC <- function(NbSeq, T, tau, Emp, pp0){
+    scan=c()
+    pvalue=c()
+    index_scan=c()
+
+    for (i in 1:NbSeq){
+        ppi=pp0[[i]]
+        result=PValue(Emp,ppi,T,tau)
+        scan=c(scan,result[1])
+        pvalue=c(pvalue,result[2])
+        index_scan=c(index_scan,result[3])
+    }
+
+    ScS_H0=data.frame(num=(1:NbSeq), scan_stat=scan, pvalue_scan=pvalue,class=c(pvalue<0.05)*1)
+    return(ScS_H0)
+}
+```
+
+```{r}
+ComputeE <- function(lambda0, lambda1){
+    E = 1
+    maxXk = floor(E*(log(lambda1/lambda0)))
+    while (maxXk < 3) {
+        E = E+1
+        maxXk = floor(E*(log(lambda1/lambda0)))
+    }
+
+    return (E)
+    }
+```
+
+```{r}
+ScoreDistribEmpiric <- function(lambda0, lambda1, n_sample, T){
+    E = ComputeE(lambda0, lambda1)
+    Score = c()
+    
+    for (i in 1:n_sample){
+        ppH0 = PoissonProcess(lambda0,T)
+        n1 = length(ppH0)
+        tbe0 = ppH0[2:n1]-ppH0[1:n1-1]
+        X = floor(E*(log(lambda1/lambda0)+(lambda0-lambda1)*tbe0))
+        Score=c(Score,X)
+    }
+    min_X = min(Score)
+    max_X = max(Score)
+
+    P_X = table(factor(Score, levels = min_X:max_X))/sum(table(Score))
+    df = data.frame("Score_X" = min(Score):max(Score), "P_X" = P_X)
+    df <- df[,-2]
+
+    return (df)
+    }
+```
+
+```{r}
+ScoreDistribElisa <- function(lambda0, lambda1, T){
+    E = ComputeE(lambda0, lambda1)
+
+    score_max = floor(E*log(lambda1/lambda0))
+
+    ## score_min compute
+    score_min_c = floor(E*log(lambda1/lambda0)+E*(lambda0-lambda1)*T)
+
+    l = seq(score_min_c,score_max,1)
+    borne_inf = (l-E*log(lambda1/lambda0))/(E*(lambda0-lambda1))
+    borne_sup = (l+1-E*log(lambda1/lambda0))/(E*(lambda0-lambda1))
+    proba.l = pexp(rate=lambda0,borne_inf)-pexp(rate=lambda0,borne_sup)
+    S = sum(proba.l)
+    new.proba.s = proba.l/S
+    df = data.frame("Score_X" = l, "P_X" = new.proba.s)
+
+    return (df)
+}
+```
+
+```{r}
+LocalScoreMC <- function(lambda0, lambda1, NbSeq, T, X_seq, P_X, tbe0){
+    E = ComputeE(lambda0, lambda1)
+    pvalue = c()
+    X = c()
+    min_X = min(X_seq)
+    max_X = max(X_seq)
+    NbSeq.NonNulles = 0
+    for (i in 1:NbSeq){
+        x = floor(E*log(dexp(tbe0[[i]], rate = lambda1)/dexp(tbe0[[i]], rate = lambda0)))
+        if (length(x)!=0){
+            X = c(X,x)
+            LS = localScoreC(x)$localScore[1]
+            daudin_result = daudin(localScore = LS, score_probabilities = P_X, sequence_length = length(x), sequence_min = min_X, sequence_max = max_X)
+            options(warn = -1) # Disable warnings print
+            pvalue = c(pvalue, daudin_result)
+            NbSeq.NonNulles = NbSeq.NonNulles + 1
+
+        }
+  }
+  LS_H0=data.frame(num=1:NbSeq.NonNulles, pvalue_scan=pvalue, class=(pvalue<0.05)*1)
+  return(LS_H0)
+}
+```
+
+```{r}
+Nb = 10
+tbe0 = vector("list",length = Nb)
+pp0 =  vector("list", length = Nb)
+pp1 =  vector("list", length = Nb)
+tbe1 = vector("list", length =  Nb)
+
+for (i in (1:Nb)) {
+    #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
+}
+
+Score = ScoreDistribEmpiric(lambda0, lambda1, Nb, T)
+LocalScoreMC(2,3,Nb,10,Score$Score_X, Score$P_X, tbe0)
+LocalScoreMC(2,3,Nb,10,Score$Score_X, Score$P_X, tbe1)
+```
+
+```{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)
+```
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+