#############################################################
# poisson.txt                                               #
# Alors voici notre premier contact avec la loi de Poisson. #
# Orsay - 20 Nov 2008                                       #
#############################################################
# Dans une expérience faite sur des rats, l’étalement sur la peau d’une crème de soin du commerce 
# peut provoquer une rougeur dans 20 cas sur 1000.
# Soit X la variable aléatoire 
# X="Nombre de rats présentant une rougeur après dépôt du produit sur N individus"
# -Q1/ Calcul d'une probabilité
#      ------------------------
# A/	Quelle loi de probabilité suit X ?
# B/	Quelle est la probabilité d'observer moins de 3 cas de réaction sur 100 étalements ? 
# 	En déduire la probabilité d'observer au moins 3 cas de réaction sur 100 étalements.
# help("dpois")
#

N=100
pelem=0.02
lambda1=(N*pelem+N*pelem*(1-pelem))/2
lambda1
etendue1=0:2
proba1=sum(dpois(etendue1,lambda1))
proba1

# pour comparer avec la distribution binomiale native :

proba2=sum(dbinom(etendue1,N,pelem))
biais=proba1-proba2
biais

#
# -Q2/ Tracer l'histogramme de la ditribution 
#      --------------------------------------

# etendue2: nombre de valeurs de la VA X considérées pour afficher le graphe
etendue2=0:40
# distribution poissonienne 
distrib1=dpois(etendue2,lambda1)
# passer le titre du graphe en variable grace à la fonction paste qui convertit et/ou concatène les éléments 
# ainsi, des paramètres numériques peuvent entrer dans la composition du titre (comme lambda)
titre=paste("X VA suivant une loi de poisson de paramètre",lambda1)
# Le résultat est un graphe dépendant des paramètre N et pelem 
#
# Recherche du maximum en Y
maxi=max(distrib1)
maxi
# D'autre part ylim=c(0,maxi), fixe les limites de l'axe y
#
plot(distrib1~etendue2,main=titre,type="h",lwd=4,col="red",xlab="X=k",ylab="P(X=k)",las=1,ylim=c(0,maxi))
# autoriser le dessin d'une nouvelle distribution sur le même graphe
par(new=TRUE)
# dessin, pour comparaison, du graphe de la binomiale native
distrib2=dbinom(etendue2,100,pelem)
plot(distrib2~etendue2,type="h",lwd=1,col="blue",xlab="X=k",ylab="P(X=k)",las=1,ylim=c(0,maxi))
# 
# Vous pouvez essayer maintenant différents couples (n, pelem) en changeant ... 
# juste les 2 premières lignes du traitement!
#####################
# FIN DU TRAITEMENT #
#####################