From 799dbef0789dfeee0c22e81b6fb6232cd737a185 Mon Sep 17 00:00:00 2001 From: Paul-Corbalan <58653590+Paul-Corbalan@users.noreply.github.com> Date: Tue, 10 May 2022 10:05:55 +0200 Subject: [PATCH 1/7] Integration scanstat.Rmd code --- Comparaison_of_methods.rmd | 141 ++++++++++++++++++++++++++++++++++++- 1 file changed, 140 insertions(+), 1 deletion(-) diff --git a/Comparaison_of_methods.rmd b/Comparaison_of_methods.rmd index 23f7a09..f529687 100644 --- a/Comparaison_of_methods.rmd +++ b/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) +``` \ No newline at end of file From 4673fb8fa7fc0445f0151c6c5ccd0a8d5e42c663 Mon Sep 17 00:00:00 2001 From: Paul-Corbalan <58653590+Paul-Corbalan@users.noreply.github.com> Date: Tue, 10 May 2022 11:06:27 +0200 Subject: [PATCH 2/7] CompareMethods used for experiences --- Comparaison_of_methods.rmd | 276 +++------------------ scanstat.Rmd | 478 ------------------------------------- 2 files changed, 34 insertions(+), 720 deletions(-) delete mode 100644 scanstat.Rmd diff --git a/Comparaison_of_methods.rmd b/Comparaison_of_methods.rmd index f529687..8ccab51 100644 --- a/Comparaison_of_methods.rmd +++ b/Comparaison_of_methods.rmd @@ -378,190 +378,6 @@ LocalScoreMC <- function(lambda0, lambda1, NbSeq, T, X_seq, P_X, tbe0){ ## 4. Experience plan for comparaison ```{r} -NbSeq = 10**2 -T = 10 - -list_of_lambda = list() -list_of_lambda[[1]] = c(1, 3) -list_of_lambda[[2]] = c(1, 4) -list_of_lambda[[3]] = c(1, 5) -list_of_lambda[[4]] = c(2, 4) -list_of_lambda[[5]] = c(2, 5) -list_of_lambda[[6]] = c(2, 6) -list_of_lambda[[7]] = c(4, 5) -list_of_lambda[[8]] = c(4, 8) -list_of_lambda[[9]] = c(4, 10) - -for (Lambda in list_of_lambda){ - lambda0 = Lambda[1] - lambda1 = Lambda[2] - Sensitivity = c() - Specificity = c() - accepted_lambda = c() - if (lambda0 < lambda1){ - - accepted_lambda = c(accepted_lambda,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) - - theoretical_results = c(rep(0,NbSeq), rep(1,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) - Emp = EmpDistrib(lambda0,n_sample,T,tau) - - X_seq = Score$Score_X - P_X = Score$P_X - - LS_H0 = LocalScoreMC(lambda0, lambda1, NbSeq, T, X_seq, P_X, tbe0) - options(warn = -1) # Disable warnings print - SS_H0 = ScanStatMC(NbSeq, T, tau, Emp, pp0) - SS_H1 = ScanStatMC(NbSeq, T, tau, Emp, pp1) - - SS_expected = c(SS_H0$class, SS_H1$class) - - #cat("Local Score:\n") - #print(summary(LS_H0)) - #cat("Scan Statistics:\n") - #print(summary(SS_H0)) - #cat("Confusion Matrix:\n") - #print(confusionMatrix(factor(LS_H0$class), factor(SS_H0$class))) - - #cat("- Elisa version:\n") - Score = ScoreDistribTheo(lambda0, lambda1, T) - Emp = EmpDistrib(lambda0,n_sample,T,tau) - - X_seq = Score$Score_X - P_X = Score$P_X - - LS_H0 = LocalScoreMC(lambda0, lambda1, NbSeq, T, X_seq, P_X, tbe0) - options(warn = -1) # Disable warnings print - - SS_H0 = ScanStatMC(NbSeq, T, tau, Emp, pp0) - - #cat("Local Score:\n") - #print(summary(LS_H0)) - #cat("Scan Statistics:\n") - #print(summary(SS_H0)) - #cat("Confusion Matrix:\n") - print(confusionMatrix(factor(theoretical_results), factor(SS_expected))) - #Sensitivity = c(Sensitivity,confusionMatrix(factor(theoretical_results), factor(SS_expected))$byClass[1]) - #Specificity = c(Specificity,confusionMatrix(factor(theoretical_results), factor(SS_expected))$byClass[2]) - - cat("---\n") - - } - titleSens=TeX(paste(r'(Sensitivity for $\lambda_0=$)', lambda0)) - plot(x=accepted_lambda,y=Sensitivity, type='l', main = titleSens) - - titleSpec=TeX(paste(r'(Specificity for $\lambda_0=$)', lambda0)) - plot(x=accepted_lambda,y=Specificity, type='l', main = titleSpec) - -} -``` - -```{r} -NbSeq = 10**2 -T = 10 -lambda0 = 2 -lambda1 = 5 -n_sample=10**4 - -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) - -theoretical_results = c(rep(0,NbSeq), rep(1,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 - } - - Emp = EmpDistrib(lambda0,n_sample,T,tau) - - SS_H0 = ScanStatMC(NbSeq, T, tau, Emp, pp0) - SS_H1 = ScanStatMC(NbSeq, T, tau, Emp, pp1) - - SS_expected = c(SS_H0$class, SS_H1$class) - - #cat("Local Score:\n") - #print(summary(LS_H0)) - #cat("Scan Statistics:\n") - #print(summary(SS_H0)) - #cat("Confusion Matrix:\n") - #print(confusionMatrix(factor(LS_H0$class), factor(SS_H0$class))) - - #cat("- Elisa version:\n") - Score = ScoreDistribTheo(lambda0, lambda1, T) - Emp = EmpDistrib(lambda0,n_sample,T,tau) - - X_seq = Score$Score_X - P_X = Score$P_X - - LS_H0 = LocalScoreMC(lambda0, lambda1, NbSeq, T, X_seq, P_X, tbe0) - options(warn = -1) # Disable warnings print - - SS_H0 = ScanStatMC(NbSeq, T, tau, Emp, pp0) - - #cat("Local Score:\n") - #print(summary(LS_H0)) - #cat("Scan Statistics:\n") - #print(summary(SS_H0)) - #cat("Confusion Matrix:\n") - print(confusionMatrix(factor(theoretical_results), factor(SS_expected))) - #Sensitivity = c(Sensitivity,confusionMatrix(factor(theoretical_results), factor(SS_expected))$byClass[1]) - #Specificity = c(Specificity,confusionMatrix(factor(theoretical_results), factor(SS_expected))$byClass[2]) - - cat("---\n") - - titleSens=TeX(paste(r'(Sensitivity for $\lambda_0=$)', lambda0)) - plot(x=accepted_lambda,y=Sensitivity, type='l', main = titleSens) - - 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){ @@ -630,66 +446,42 @@ CompareMethods <- function(lambda0, lambda1, NbSeq, T, tau){ ``` ```{r} -NbSeq = 100 +NbSeq = 10**2 T = 10 tau = 2 -lambda0 = 0.3 -lambda1 = 0.5 -result1 = CompareMethods(lambda0, lambda1, NbSeq, T, tau) +list_of_lambda = list() +list_of_lambda[[1]] = c(1, 3) +list_of_lambda[[2]] = c(1, 4) +list_of_lambda[[3]] = c(1, 5) +list_of_lambda[[4]] = c(2, 4) +list_of_lambda[[5]] = c(2, 5) +list_of_lambda[[6]] = c(2, 6) +list_of_lambda[[7]] = c(4, 5) +list_of_lambda[[8]] = c(4, 8) +list_of_lambda[[9]] = c(4, 10) -lambda0 = 0.01 -lambda1 = 1 +i = 1 +legend_list = c() -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) -``` \ No newline at end of file +for (Lambda in list_of_lambda){ + lambda0 = Lambda[1] + lambda1 = Lambda[2] + result = CompareMethods(lambda0, lambda1, NbSeq, T, tau) + title_ROC = TeX(paste(r'(ROC curve for several values of $\lambda_0$ and $\lambda_1$)')) + + perfSS = result[1] + perfLS = result[2] + + + plot(perfSS$performance.SS, lty=1, col=i, lwd = 2) + par(new=T) + plot(perfLS$performance.LS, lty=2, col=i,lwd = 2) + + legend_list=c(legend_list, paste(c("lambda0 = ", lambda0, ", lambda1 = ", lambda1), collapse = "")) + + i=i+1 +} + +legend(0.5, 0.3, legend=legend_list, col=1:length(list_of_lambda), lty=1, cex=0.9,lwd=4, box.lty=0) +``` diff --git a/scanstat.Rmd b/scanstat.Rmd deleted file mode 100644 index 1278839..0000000 --- a/scanstat.Rmd +++ /dev/null @@ -1,478 +0,0 @@ ---- -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(ppH0bisppH1.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) Date: Tue, 10 May 2022 11:21:14 +0200 Subject: [PATCH 3/7] 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) Date: Tue, 10 May 2022 11:30:40 +0200 Subject: [PATCH 4/7] Add tau_scan, tau_H1 in CompareMethods --- Comparaison_of_methods.rmd | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Comparaison_of_methods.rmd b/Comparaison_of_methods.rmd index 92e04d0..9eb4225 100644 --- a/Comparaison_of_methods.rmd +++ b/Comparaison_of_methods.rmd @@ -380,7 +380,7 @@ LocalScoreMC <- function(lambda0, lambda1, NbSeq, T, X_seq, P_X, tbe0){ ## 4. Experience plan for comparaison ```{r} -CompareMethods <- function(lambda0, lambda1, NbSeq, T, tau){ +CompareMethods <- function(lambda0, lambda1, NbSeq, T, tau_scan, tau_H1){ if (lambda0 < lambda1){ cat("For T = ", T, ", Nb = ", NbSeq, ", lambda0 = ", lambda0, " and lambda1 = ", lambda1, ":\n", sep = "") @@ -398,7 +398,7 @@ CompareMethods <- function(lambda0, lambda1, NbSeq, T, tau){ tbe0[[i]] = tbei #Simulation for sequences under H1 - ppj1 = SimulationH1(lambda0, lambda1, T, 3) + ppj1 = SimulationH1(lambda0, lambda1, T, tau_H1) nj = length(ppj1) pp1[[i]] = ppj1 tbej = ppj1[2:nj]-ppj1[1:nj-1] @@ -412,9 +412,9 @@ CompareMethods <- function(lambda0, lambda1, NbSeq, T, tau){ LS_obtained = c(LS_H0$class, LS_H1$class) options(warn = -1) - Emp = EmpDistrib(lambda0,10**4,T,tau) - SS_H0 = ScanStatMC(NbSeq, T, tau, Emp, pp0) - SS_H1 = ScanStatMC(NbSeq, T, tau, Emp, pp1) + # Emp = EmpDistrib(lambda0, 10**4, T, tau_scan) + SS_H0 = ScanStatMC(NbSeq, T, tau_scan, Emp, pp0) + SS_H1 = ScanStatMC(NbSeq, T, tau_scan, Emp, pp1) SS_obtained = c(SS_H0$class, SS_H1$class) From 0054d077be4d9c6752e6e8f9800da8bacf0f4f85 Mon Sep 17 00:00:00 2001 From: Paul-Corbalan <58653590+Paul-Corbalan@users.noreply.github.com> Date: Tue, 10 May 2022 18:45:12 +0200 Subject: [PATCH 5/7] Add NbSeq.NonNulles to LocalScoreMC --- Comparaison_of_methods.rmd | 40 +++++++++++++++++++++----------------- 1 file changed, 22 insertions(+), 18 deletions(-) diff --git a/Comparaison_of_methods.rmd b/Comparaison_of_methods.rmd index 9eb4225..0aaf67c 100644 --- a/Comparaison_of_methods.rmd +++ b/Comparaison_of_methods.rmd @@ -356,24 +356,28 @@ plot_graph_distrib_score(distrib_score_theo, distrib_score_mc) ### 3.2. Local score calculation ```{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) - - for (i in 1:NbSeq){ - x = floor(E*log(dexp(tbe0[[i]], rate = lambda1)/dexp(tbe0[[i]], rate = lambda0))) - 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) + E = ComputeE(lambda0, lambda1) + + pvalue = 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){ + 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, pvalue_scan=pvalue, class=as.numeric(pvalue<0.05)) + LS_H0=data.frame(num=1:NbSeq.NonNulles, pvalue_scan=pvalue, class=as.numeric(pvalue<0.05)) return(LS_H0) } ``` @@ -451,7 +455,7 @@ CompareMethods <- function(lambda0, lambda1, NbSeq, T, tau_scan, tau_H1){ ``` ```{r} -NbSeq = 10**2 +NbSeq = 10**4 T = 10 tau = 2 From 980d36c17511fb554e781383aef5ae324fd7720f Mon Sep 17 00:00:00 2001 From: njulia1 Date: Sat, 14 May 2022 17:49:21 +0200 Subject: [PATCH 6/7] Update Comparaison_of_methods.rmd Modification of LOCAL SCORE MC function and hypothesis test --- Comparaison_of_methods.rmd | 49 ++++++++++++++++++-------------------- 1 file changed, 23 insertions(+), 26 deletions(-) diff --git a/Comparaison_of_methods.rmd b/Comparaison_of_methods.rmd index 92e04d0..9c0fabd 100644 --- a/Comparaison_of_methods.rmd +++ b/Comparaison_of_methods.rmd @@ -356,24 +356,25 @@ plot_graph_distrib_score(distrib_score_theo, distrib_score_mc) ### 3.2. Local score calculation ```{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) - - for (i in 1:NbSeq){ - x = floor(E*log(dexp(tbe0[[i]], rate = lambda1)/dexp(tbe0[[i]], rate = lambda0))) - 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) + 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, pvalue_scan=pvalue, class=as.numeric(pvalue<0.05)) + LS_H0=data.frame(num=1:NbSeq.NonNulles, pvalue_scan=pvalue, class=(pvalue<0.05)*1) return(LS_H0) } ``` @@ -451,20 +452,16 @@ CompareMethods <- function(lambda0, lambda1, NbSeq, T, tau){ ``` ```{r} -NbSeq = 10**2 +NbSeq = 10**4 T = 10 tau = 2 list_of_lambda = list() list_of_lambda[[1]] = c(1, 3) -list_of_lambda[[2]] = c(1, 4) -list_of_lambda[[3]] = c(1, 5) -list_of_lambda[[4]] = c(2, 4) -list_of_lambda[[5]] = c(2, 5) -list_of_lambda[[6]] = c(2, 6) -list_of_lambda[[7]] = c(4, 5) -list_of_lambda[[8]] = c(4, 8) -list_of_lambda[[9]] = c(4, 10) +list_of_lambda[[2]] = c(2, 6) +list_of_lambda[[3]] = c(4, 7) +list_of_lambda[[4]] = c(2, 9) + i = 1 legend_list = c() From e3eba7b442f96d82f905ce776534956042dec3cd Mon Sep 17 00:00:00 2001 From: njulia1 Date: Sun, 15 May 2022 15:38:34 +0200 Subject: [PATCH 7/7] Update Comparaison_of_methods.rmd Computation of the error for the simulation of sequences under H_0 : change of values of lambda --- Comparaison_of_methods.rmd | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/Comparaison_of_methods.rmd b/Comparaison_of_methods.rmd index 9c0fabd..a0425d1 100644 --- a/Comparaison_of_methods.rmd +++ b/Comparaison_of_methods.rmd @@ -142,8 +142,9 @@ NbSeqH0=10000 NbSeqH1=NbSeqH0 DataH0=vector("list") DataH1=vector("list") -lambda0=2 +lambda0=1 lambda1=5 + T=10 tau=1 @@ -301,9 +302,6 @@ ScoreDistribTheo <- 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) @@ -334,10 +332,6 @@ distrib_score_mc = ScoreDistribEmpiric(2,3,10000,T) distrib_score_theo = ScoreDistribTheo(2,3,T) plot_graph_distrib_score <- function(distrib_score_theo, distrib_score_mc){ - # length(distrib_score_mc[,2]) - # length(distrib_score_theo[,2]) - - #diff_distrib_score=abs(distrib_score_mc[,2]-distrib_score_theo[,2]) #par(mfrow = c(1,2)) @@ -407,13 +401,13 @@ CompareMethods <- function(lambda0, lambda1, NbSeq, T, tau){ } #cat("- Empiric version:\n") - Score = ScoreDistribEmpiric(lambda0, lambda1, 10**4, T) + Score = ScoreDistribEmpiric(lambda0, lambda1, 10**5, 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,10**4,T,tau) + Emp = EmpDistrib(lambda0,10**5,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)