From fd3ab9ac2eb4ed92ba4dbfda25c78ff6419dd02f Mon Sep 17 00:00:00 2001 From: elisaduz Date: Sat, 19 Mar 2022 11:56:36 +0100 Subject: [PATCH] Update Dataset_study.rmd Ajout_Code_Troncage --- Dataset_study.rmd | 77 ++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 76 insertions(+), 1 deletion(-) diff --git a/Dataset_study.rmd b/Dataset_study.rmd index d26ea42..bfafa18 100644 --- a/Dataset_study.rmd +++ b/Dataset_study.rmd @@ -279,7 +279,7 @@ library(stringr) E=10 lambda0=1 lambda1=2 - +max=0 while (exp(-lambda0*max) > 10^(-9)){ max = max + 1 @@ -292,6 +292,11 @@ t=seq(0,max,0.1) x=floor(E*log(lambda1/lambda0)+(lambda0-lambda1)*t) head(x) range(x) + + + + + s=seq(-15,6,0.1) ##tronquer la queue des x (x négatifs) A=1/(lambda0-lambda1) @@ -328,6 +333,76 @@ p ``` +```{r} + +library(stringr) + +E=10 +lambda0=1 +lambda1=2 +max=0 + +while (exp(-lambda0*max) > 10^(-9)){ + max = max + 1 +} +max + + # en fonction de lambda0 tq P(Expo(lambda +#0)>=max)<10^-9 +t=seq(0,max,0.1) +probaseuilmin = 10^-9 +x=floor(E*log(lambda1/lambda0)+(lambda0-lambda1)*t) +head(x) +range(x) + + +s=seq(-15,6,0.1) +##tronquer la queue des x (x négatifs) +A=1/(lambda0-lambda1) +B=A*log(lambda1/lambda0) +proba.s=lambda0*exp(-lambda0*(A*s-B)) +proba.s = proba.s/sum(proba.s) +barplot(proba.s) + +# Troncage à un score minimal +minXk = as.numeric(x[max(which(proba.s<(probaseuilmin)))]) # On définit la classe pour proba < probaseuilmin +names(proba.s) = as.character(x) +proba.s[as.character(minXk)] = sum(proba.s[which(proba.s=minXk)] # On ne garde que les scores supérieurs à minXk +subtitle = paste("lambda0=",lambda0,";lambda1=",lambda1,";E=",E,sep="") +barplot(proba.s, col="steelblue",xlab="Score",ylab="Probabilité", main = paste("Probabilité d'apparition de chaque score\nLoi géométrique : ",subtitle,sep="")) + +barplot(proba.s) +##trouver 1/lambda en fonction intervalle seq +## ecart ou rapport entre les lambda qui est le + important ? + +for (lambda0 in c(5,50,500,5000,50000) ) { + for (lambda1 in c(2,20,200,2000,20000) ) { + probaS_vect=vector(length=length(s)) + + for (j in 1:length(x)) + { + + + A=1/(lambda0-lambda1) + B=A*log(lambda1/lambda0) + + probaS=floor(E*lambda0*exp(-lambda0*(A*x[j]-B))) + # print(s) + probaS_vect[j]=probaS + } + +probaS_vect = probaS_vect/sum(probaS_vect) +print(probaS_vect) +p <- barplot(probaS_vect, main=str_c("lambda0=", as.character(lambda0), "et", "lambda1=", as.character(lambda1), sep=" " )) +p +} +} + + + +``` + Import data of rainfall in France every 3 hours.