Newer
Older
# objective:
# author : I.Sanchez
# creation: 14/01/2021
#---------------------------------------------------------
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
#----------------------------------------------------------------
# Fonctions de l'algorithme de détection de rupture
#----------------------------------------------------------------
#--------------------------------------------------------------------
# objective: quelques fonctions pour detection ruptures
# author : I.Sanchez
# creation: 16/01/2023
# update : 23/02/2023
#-------------------------------------------------------------------
# sim_changepointPP: simulation de rupture
# loglikelihood: calcul vraisemblance
# reconstructionlambdamu: reconstruction de lambda et mu sachant S
# reconstructionS_avecL: reconstruction de S sachant lambda et mu
# et avec calcul de vraisemlbance
# reconstructionS_avecBH: reconstruction de S sachant lambda et mu et
# calcul p.value corrigée par BH
#-------------------------------------------------------------------
#' Fonction de simulation de ruptures suivant processus de Poisson
#' fonction qui simule le processus de poisson composee avec une rupture
#'
#' @param tau0 numeric
#' @param lambda numeric, temps moyen entre 2 chevauchements
#' @param mu numeric, temps moyen de chevauchements
#' @param t_final numeric, temps final de l'experimentation
#'
#' @return
#' @export
#'
#' @examples
sim_changepointPP <- function(tau0, lambda, mu, t_final){
#fonction qui simule le processus de poisson composee avec une rupture
N_minus <- rpois(n= 1, lambda = lambda[1]*tau0)
N_plus <- rpois(n= 1, lambda = lambda[2]*(t_final-tau0))
time_minus <- sort(runif(n = N_minus, min =0, max = tau0))
time_plus <- sort(runif(n = N_plus, min =tau0, max = t_final))
marks_minus <- rexp(N_minus, rate = mu[1])
marks_plus <- rexp(N_plus, rate = mu[2])
return(tibble(time = c(time_minus, time_plus),
marks = c(marks_minus, marks_plus),
state = as.factor(c(rep(0, N_minus), rep(1,N_plus)))
)
)
}
#-------------------------------------------------------------------
#' Calcul la vraisemblance pour un individu
#'
#' @param lambda numeric, un vecteur de taille 2 - tps moyen entre 2 chev
#' @param mu numeric, un vecteur de 2 elements - tps moyen des chev
#' @param S numeric, temps de rupture
#' @param k num, le nombre de sauts avant S
#' @param ntot num, le nombre total de sauts
#' @param U num, le vecteur des taille de sauts
#'
#' @return
#' @export
#'
#' @examples
loglikelihood_tmp <-function(lambda,mu,S,k,ntot, U,t_final){
# section 1. calcul de la vraisemblance doc XXX.pdf
# U: durée de chevauchement
# S: rupture
# k: nombre de chev avant la rupture
# ntot: nombre total de chevauchements
if (k==0){
L = ntot * log(lambda[2]*mu[2]) + (lambda[2] - lambda[1])*S
- t_final*lambda[2] - mu[2] * sum(U,na.rm=TRUE)
} else if (k==ntot){
L = ntot * log(lambda[1]*mu[1]) + (lambda[2] - lambda[1])*S
- t_final*lambda[2] - mu[1] * sum(U[1:k],na.rm=TRUE)
} else if (k != 0 & k != ntot){
L = k*log(lambda[1]*mu[1] / (lambda[2]*mu[2]) ) +
ntot * log(lambda[2]*mu[2]) + (lambda[2] - lambda[1])*S
- t_final*lambda[2] - mu[1] * sum(U[1:k],na.rm=TRUE)
- mu[2] * sum(U[(k+1):ntot],na.rm=TRUE)
}
return(L)
}
loglikelihood <-function(lambda,mu,S,k,ntot, U,t_final){
max(loglikelihood_tmp(lambda,mu,S,k,ntot, U,t_final),
loglikelihood_tmp(lambda,mu,S,k-1,ntot, U,t_final))
} else {
(lambda[2] - lambda[1])*S- t_final*lambda[2]
}
}
#-----------------------------------------------------------
#' Title
#'
#' @param datain input data.frame
#' @param S numeric, rupture
#'
#' @return
#' @export
#'
#' @examples
reconstructionlambdamu <-function(datain,S,t_final){
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
# Initialisation
nb_indiv=datain$id[length(datain$id)]
SommeS=SommeSautAvant=SommeSautApres=SommeUAvant=SommeUApres=SommeSaut=0
for (ind in 1:nb_indiv){
# ntot: nombre total de chevauchements => nombre de lignes pr 1 ind
ewe_ind=datain[datain$id==ind,]
ntot=dim(ewe_ind)[1]
SommeSaut=SommeSaut+ntot
if (ntot>0){
sauts=ewe_ind$time
n=sum(sauts < S[ind])
SommeSautAvant=SommeSautAvant + n
SommeSautApres=SommeSautApres + sum(sauts >= S[ind],na.rm=TRUE)
SommeS=SommeS + S[ind]
U =ewe_ind$marks
if(n==0){
SommeUApres=SommeUApres + sum(U[1:ntot],na.rm=TRUE)
} else if(n==ntot){
SommeUAvant=SommeUAvant + sum(U[1:ntot],na.rm=TRUE)
} else {
SommeUAvant=SommeUAvant + sum(U[1:n],na.rm=TRUE)
SommeUApres=SommeUApres + sum(U[(n+1):ntot],na.rm=TRUE)
}
}
}
# lambdan[1]: tps moyen entre 2 chev si pas en chaleur (avant rupture)
# lambdan[2]: tps moyen entre 2 chev si en chaleur (après rupture)
# mun[1]: tps moyen des chev si pas en chaleur (avant rupture)
# mun[2]: tps moyen des chev si en chaleur (après rupture)
if (SommeSautAvant==0){SommeUAvant=1}
if (SommeSautApres==0){SommeUApres=1} #garde fou pour dire 0/0=0
if (max(SommeSautAvant/SommeS,0.01) < max(SommeSautApres/SommeS,0.01)){
# mes deux max sont des gardes fous pour ne pas avoir lambda=0
return(c(max(SommeSautAvant/SommeS,0.01),
max(SommeSautApres/(nb_indiv*t_final - SommeS),0.01),
max(SommeSautAvant/SommeUAvant,0.01),
max(SommeSautApres/SommeUApres,0.01))
)
} else {
return(c(max(SommeSaut/(nb_indiv*t_final),0.01),
max(SommeSaut/(nb_indiv*t_final)/SommeS,0.01),
max(SommeSautAvant/SommeUAvant,0.01),
max(SommeSautApres/SommeUApres,0.01))
)
}
}
#-----------------------------------------------------------
#' reconstruction des temps de Sauts avec la vraisemblance
#'
#' @param lambdan numeric, un vecteur de taille 2 - tps moyen entre 2 chev
#' @param mun numeric, un vecteur de 2 elements - tps moyen des chev
#' @param ewe dataset pour le MC
#'
#' @details S tps de sauts (rupture) et U duree du chevauchement
#' @return
#' @export
#'
#' @examples
reconstructionS_avecL <- function(lambdan,mun,ewe,t_final){
S_reconstruit=c()
l0=l1=c()
nb_indiv=ewe$id[length(ewe$id)]
for (ind in 1:nb_indiv){
ewe_ind=ewe[ewe$id==ind,]
# nombre total de chevauchements
ntot=dim(ewe_ind)[1]
if (ntot==0){
S_reconstruit=c(S_reconstruit,t_final)
l0=c(l0,1)
l1=c(l1,1) # le 1 ici on s'en fout puisque l1=l0 et on regardera la diff
} else {
ntot=dim(ewe_ind)[1]
# Initialisation avec k un saut, k=1 <=> saut 1
k=1
kopt=k
S = ewe_ind$time[k]
U = ewe_ind$marks
# Vraisemblance optimale
likelihoodopt<-loglikelihood(lambdan,mun,S,k,ntot, U,t_final)
Sopt=S
# Optimisation si ntot > 1
if (ntot > 1){
for (k in 2:ntot){
S=ewe_ind$time[k]
likelihoodk<-loglikelihood(lambdan,mun,S,k,ntot, U,t_final)
if (likelihoodk > likelihoodopt){
likelihoodopt=likelihoodk
kopt=k
Sopt=S
}
}
# BC: ici il faudrait ajouter une ligne pour ajouter t_final
# mais cela me semble inutile avec la couche BH...
}
S_reconstruit=c(S_reconstruit,Sopt)
l0=c(l0,loglikelihood(lambdan,mun,t_final,ntot,ntot, U,t_final))
l1=c(l1,likelihoodopt)
}
}
# retourne une liste avec les ruptures reconstruites
# et la différence des vraisemblances l1 - l0
return(list(S_reconstruit,l1-l0))
}
#-----------------------------------------------------------
#' Création d'ech par bootstrap paramétrique pour calculer les pvalues
#' avec correction de la multiplicité (FDR=BH) par MC
#'
#' @param datain
#' @param alpha
#' @param nb_MC
#' @param lambdan
#' @param mun
#'
#' @return
#' @export
#'
#' @examples
reconstructionS_avecBH<-function(datain,alpha,nb_MC,lambdan,mun,t_final){
nb_indiv=nb_MC
list_individuals <- lapply(1:nb_indiv, function(i){
tau0 = t_final
ewe_MC <- sim_changepointPP(tau0 = tau0 , lambda = lambdan,
mu = mun, t_final = t_final)
ewe_MC <- ewe_MC %>% mutate(id = i, tau = tau0)
})
ewe_MC <- do.call('rbind', list_individuals)
Sl_data=reconstructionS_avecL(lambdan,mun,datain,t_final)
Sl_MC=reconstructionS_avecL(lambdan,mun,ewe_MC,t_final)
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
l_data=Sl_data[[2]]
l_MC=Sl_MC[[2]]
S=Sl_data[[1]]
nb_indiv=datain$id[length(datain$id)]
pvalue=c()
for (ind in 1:nb_indiv){
#cat("reconstructionS_avecBH - animal:",ind,"\n")
# calcul des pvalues
pvalue=c(pvalue, sum(l_MC > l_data[ind]) / nb_MC)
}
pvalueclasse=sort(pvalue, index.return = TRUE)
# Correction multiplicité par FDR: méthode Benjamini-Hochberg
#--------------------------------------------------------------
# alpha == q dans FDR/pdf section 1
critBH<-numeric()
BH = 0 # passe à 1 dès que j'arrête d'attribuer des ruptures
for (ind in 1:nb_indiv){
if (BH==0 & (pvalueclasse$x[ind] > ind*alpha/nb_indiv)){
BH=1
}
if (BH==1){
S[pvalueclasse$ix[ind]]=t_final
}
critBH[ind]<-BH
}
# liste: S (rupture), les p.values, l'ordre des p.values et le critère BH
return(list(S,pvalue,pvalueclasse,critBH))
}
#----------------------------------------------------------------
# Appel des fonctions de l'algorithme de détection de ruptures
#----------------------------------------------------------------
#' @param datain input data.frame
#' @param lambda0 vecteur numérique de taille 2, temps entre 2 chevauchements
#' @param mu0 vectuer numérique de taille 2, durée des chevauchements
#' @param nprec numeric, n total de convergence vers l'optimum (n precision)
#' @param nb_MC numeric, nombre d'echantillons MC simulés
#'
#' @return
#' @export
#'
#' @examples
detectS_BH<- function(datain,lambda0,mu0,nprec,alpha,nb_MC,t_final){
# boucle sur chaque brebis, algorithme de détection de rupture
# avec correction multiplicité par méthode BH
#' @param nprec numeric, n total de convergence vers l'optimum (n precision)
#' @param nb_MC numeric, nombre d'echantillons MC simulés
# Initialisation
lambdan<-lambda0
mun<-mu0
#nprec<-100
#alpha<-0.05
#nb_MC<-100
for (i in 1:nprec){
# Boucle de convergence vers l'optimum
outBH=reconstructionS_avecBH(datain,alpha,nb_MC,lambdan,mun,t_final)
names(outBH)<-c("S","p.value","p.value.classe","BH")
S=outBH[[1]]
Param=reconstructionlambdamu(datain,S,t_final)
lambdan=Param[1:2]
mun=Param[3:4]
}
# output list: outBH (S, pval et pvalclasse) et lambdan et mun
output<-list(outBH,lambdan,mun)
names(output)<-c("outBH","lambda_n","mu_n")
return(output)
}