CompareMethods used for experiences
This commit is contained in:
Paul-Corbalan
parent
799dbef078
commit
4673fb8fa7
|
|
@ -378,190 +378,6 @@ LocalScoreMC <- function(lambda0, lambda1, NbSeq, T, X_seq, P_X, tbe0){
|
||||||
|
|
||||||
## 4. Experience plan for comparaison
|
## 4. Experience plan for comparaison
|
||||||
```{r}
|
```{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){
|
CompareMethods <- function(lambda0, lambda1, NbSeq, T, tau){
|
||||||
if (lambda0 < lambda1){
|
if (lambda0 < lambda1){
|
||||||
|
|
||||||
|
|
@ -630,66 +446,42 @@ CompareMethods <- function(lambda0, lambda1, NbSeq, T, tau){
|
||||||
```
|
```
|
||||||
|
|
||||||
```{r}
|
```{r}
|
||||||
NbSeq = 100
|
NbSeq = 10**2
|
||||||
T = 10
|
T = 10
|
||||||
tau = 2
|
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
|
i = 1
|
||||||
lambda1 = 1
|
legend_list = c()
|
||||||
|
|
||||||
result2 = CompareMethods(lambda0, lambda1, NbSeq, T, tau)
|
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$)'))
|
||||||
|
|
||||||
lambda0 = 1
|
perfSS = result[1]
|
||||||
lambda1 = 1.1
|
perfLS = result[2]
|
||||||
|
|
||||||
result3 = CompareMethods(lambda0, lambda1, NbSeq, T, tau)
|
|
||||||
|
|
||||||
|
|
||||||
lambda0 = 0.9
|
plot(perfSS$performance.SS, lty=1, col=i, lwd = 2)
|
||||||
lambda1 = 2
|
par(new=T)
|
||||||
|
plot(perfLS$performance.LS, lty=2, col=i,lwd = 2)
|
||||||
|
|
||||||
result4 = CompareMethods(lambda0, lambda1, NbSeq, T, tau)
|
legend_list=c(legend_list, paste(c("lambda0 = ", lambda0, ", lambda1 = ", lambda1), collapse = ""))
|
||||||
|
|
||||||
|
i=i+1
|
||||||
|
}
|
||||||
|
|
||||||
title_ROC = TeX(paste(r'(ROC curve for several values of $\lambda_0$ and $\lambda_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)
|
||||||
|
|
||||||
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)
|
|
||||||
```
|
```
|
||||||
478
scanstat.Rmd
478
scanstat.Rmd
|
|
@ -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(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)
|
|
||||||
```
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
Loading…
Reference in New Issue