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