paul-corbalan
/
Scan-Statistics-Project-4Y-INSA
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Integration of tau in compare and PValue debug

This commit is contained in:
Paul-Corbalan 2022-05-10 11:21:14 +02:00
parent 4673fb8fa7
commit 864fd506a1
1 changed files with 22 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

39
Comparaison_of_methods.rmd
View File

@ -12,6 +12,7 @@ library("latex2exp")
library("Rcpp") library("Rcpp")
library("caret") library("caret")
library("ROCR") library("ROCR")
library("pROC")
``` ```
## 1. Proposition for simulations under $\mathcal{H}_1$ ## 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`: Compute $p$-value for scan statistic of `ppH1` with `Emp`:
```{r} ```{r}
PValue <- function(Emp,ppH, T, tau){ PValue <- function(Emp,ppH, T, tau){
scanH1=ScanStat(ppH,T,tau)[2] SS = ScanStat(ppH,T,tau)
index_scanH1=ScanStat(ppH,T,tau)[1] scanH = SS[2]
index=Emp$index_scan index_scanH = SS[1]
index = Emp$index_scan
n=length(index) n=length(index)
if (scanH1< min(Emp$index_scan)){ if (scanH< min(Emp$index_scan)){
return (c(scanH1,1,index_scanH1)) return (c(scanH,1,index_scanH))
} else{ } else{
if(min(Emp$index_scan)<scanH1 && scanH1<=max(Emp$index_scan)){ if(min(Emp$index_scan)<scanH && scanH<=max(Emp$index_scan)){
return(c(scanH1,1-Emp$cdf[scanH1-min(Emp$index_scan)+1],index_scanH1)) return(c(scanH,1-Emp$cdf[scanH-min(Emp$index_scan)],index_scanH))
} else{return (c(scanH1,0,index_scanH1))}} } else{return (c(scanH,0,index_scanH))}}
} }
``` ```
### 2.2. Simulation under $\mathcal{H}_0$ and computation of p-values ### 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){ if (lambda0 < lambda1){
cat("For T = ", T, ", Nb = ", NbSeq, ", lambda0 = ", lambda0, " and lambda1 = ", lambda1, ":\n", sep = "") 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) pp0 = vector("list", length = NbSeq)
pp1 = vector("list", length = NbSeq) pp1 = vector("list", length = NbSeq)
tbe1 = vector("list", length = NbSeq) tbe1 = vector("list", length = NbSeq)
@ -396,7 +398,7 @@ CompareMethods <- function(lambda0, lambda1, NbSeq, T, tau){
tbe0[[i]] = tbei tbe0[[i]] = tbei
#Simulation for sequences under H1 #Simulation for sequences under H1
ppj1 = SimulationH1(lambda0, lambda1, T, tau) ppj1 = SimulationH1(lambda0, lambda1, T, 3)
nj = length(ppj1) nj = length(ppj1)
pp1[[i]] = ppj1 pp1[[i]] = ppj1
tbej = ppj1[2:nj]-ppj1[1:nj-1] tbej = ppj1[2:nj]-ppj1[1:nj-1]
@ -404,13 +406,13 @@ CompareMethods <- function(lambda0, lambda1, NbSeq, T, tau){
} }
#cat("- Empiric version:\n") #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_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_H1 = LocalScoreMC(lambda0, lambda1, NbSeq, T, Score$Score_X, Score$P_X, tbe1)
LS_obtained = c(LS_H0$class, LS_H1$class) LS_obtained = c(LS_H0$class, LS_H1$class)
options(warn = -1) 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_H0 = ScanStatMC(NbSeq, T, tau, Emp, pp0)
SS_H1 = ScanStatMC(NbSeq, T, tau, Emp, pp1) SS_H1 = ScanStatMC(NbSeq, T, tau, Emp, pp1)
SS_obtained = c(SS_H0$class, SS_H1$class) 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))) 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), print(confusionMatrix(as.factor(SS_obtained), as.factor(theoretical_results_SS),
dnn = c("Prediction", "Reference"))$table) 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") cat("--- Confusion matrix for local score method --- \n")
theoretical_results_LS = c(rep(0,length(LS_H0$num)), rep(1,length(LS_H1$num))) 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), print(confusionMatrix(as.factor(LS_obtained), as.factor(theoretical_results_LS),
dnn = c("Prediction", "Reference"))$table) dnn = c("Prediction", "Reference"))$table)
#cat("--- Coube ROC associé")
title_ROC = TeX(paste(r'(ROC curve for $H_0: \lambda_0=$)', lambda0, title_ROC = TeX(paste(r'(ROC curve for $H_0: \lambda_0=$)', lambda0,
r'(against $H_1: \lambda_0=$)', lambda1)) r'(against $H_1: \lambda_0=$)', lambda1))
pred.SS = prediction(theoretical_results_SS,SS_obtained) pred.SS = prediction(theoretical_results_SS,SS_obtained)
pred.LS = prediction(theoretical_results_LS,LS_obtained) pred.LS = prediction(theoretical_results_LS,LS_obtained)
perf.SS = performance(pred.SS,"tpr", "fpr") perf.SS = performance(pred.SS,"tpr", "fpr")
perf.LS = performance(pred.LS,"tpr", "fpr") perf.LS = performance(pred.LS,"tpr", "fpr")
#plot(perf.SS, lty=1, col="coral")
par(new=T) 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") cat("-----------------------------------\n")
options(warn = -1)
result <- c('performance.SS'= perf.SS,'performance.LS'= perf.LS) result <- c('performance.SS'= perf.SS,'performance.LS'= perf.LS)
return(result) return(result)