Exercice R et Stat Des
Christophe Ambroise
2024-03-04
Simulation de lois
Loi uniforme: simulation
Simuler 100 réalisation d’une loi uniforme
set.seed(2)
x<-runif(100,5,10)
plot(x)
et compter combien de valeurs simulées sont au dessus d’un seuil
compteur<-0
for(value in x) if(value>=7) compteur<-compteur+1
print(compteur)## [1] 56En version courte
# Version courte
sum(x>=7)## [1] 56Loi Gaussienne
- Simuler 300 valeurs de QI qui suivent une loi normale \[ QI \sim \mathcal{N} (\mu=100, \sigma=15) \]
- Estimer la proportion de valeurs
- supérieures à 130
- inférieures à 100
QI<-rnorm(300,mean=100,sd=15)
mean(QI<100)## [1] 0.4666667mean(QI>130)## [1] 0.03333333- Calculer la valeur théorique de ces proportions
\[ \begin{align} P(QI>130) & = P(\frac{QI-100}{15}>\frac{130-100}{15})\\ & = P(\frac{QI-100}{15}>2)\\ & = 1 - P(\frac{QI-100}{15}\leq 2)\\ & = 1 -P(U \leq 2)\\ & = 1- \phi(2)\\ & = \phi(-2) \end{align} \] avec \(U \sim \mathcal N(0,1)\).
print(pnorm(130,mean=100, sd=15,lower.tail = FALSE))## [1] 0.02275013print(1 - pnorm(2))## [1] 0.02275013print(pnorm(-2))## [1] 0.02275013- Utiliser ces calculs pour illustrer la loi des grands nombres
La loi des grands nombres nous dit que
\[ \frac{1}{n}\sum_{j=1}^n \mathbb{I}_{QI_j >130} \stackrel{p.s.}{\longrightarrow} E(\mathbb I_{QI_j}>130) = p_{130} \] Si on augmente la taille d’échantillon et que l’on calcule la proportion empirique:
taille.echantillon<-c(10,100,1000,10000,100000)
for (n in taille.echantillon){print(mean(rnorm(n=n,mean=100,sd=15)>130))}## [1] 0.1
## [1] 0.03
## [1] 0.025
## [1] 0.0238
## [1] 0.02265Prétraitement
Outliers (données abberantes)
- Soit $k {1,2,…} et X une variable normale. Calculer \[ \forall k, \ P(X \in [\mu-k\sigma,\mu+k\sigma]) \]
# Simulation de 10,000 tailles
x<-rnorm(10000,mean=0.5,sd=2/100)
hist(x)
x<- c(x,5/100,5)Estimons la moyenne et l’écart type de l’échantillon et prenons k=4
# Simulons 1000 taille de nourissons suivant une loi normale
x<-rnorm(1000,mean=50,sd=2)
# Ajoutons deux erreurs (1 bébé trop petit, 1 trop grand)
c(x,4.8,510)->x
# Estimons la moyenne et l'écart de l'échantillon des 1002 bébés
mu.hat<-mean(x)
sd.hat<-sqrt(var(x))
k<-4
print(outliers.min<- mu.hat - k*sd.hat)## [1] -8.577834print(outliers.max<- mu.hat + k*sd.hat)## [1] 109.4141x[(x<outliers.min)|(x> outliers.max)]## [1] 510Calculons la moyenne et la variance tronquée en utilisant un paramètre trim qui donne le pourcentage d’observations extrêmes non pris en compte…
trim<-0.01
n<-length(x)
k<-round(n*trim)
x.trimmed <-sort(x)[(k+1):(n-k)]
mu.hat<-mean(x.trimmed)
sd.hat<-sqrt(var(x.trimmed))
sd.trim<-function(x,trim=0.05){
n<-length(x)
k<-round(n*trim)
x.trimmed <-sort(x)[(k+1):(n-k)]
sd.hat<-sqrt(var(x.trimmed))
}
k<-4
print(outliers.min<- mu.hat - k*sd.hat)## [1] 42.18132print(outliers.max<- mu.hat + k*sd.hat)## [1] 57.81216x[(x<outliers.min)|(x> outliers.max)]## [1] 58.3388 4.8000 510.0000which((x<outliers.min)|(x> outliers.max))## [1] 49 1001 1002Fonction
Ecrire une fonction qui prend en argument un échantillon, la variable trim qui indique le pourcentage d’observations à supprimer, et qui rende en sortie, la moyenne tronquée et l’écart type tronqué.
stat.trimmed<-function(x,trim=0.05){
n<-length(x)
k<-round(n*trim)
x.trimmed<-sort(x)[(k+1):(n-k)]
return(list(mean.trimmed=mean(x.trimmed),
sd.trimmed=sqrt(var(x.trimmed))))
}
# simulons un échantillon
x<-c(rnorm(100,50,2),600)
print(stat.trimmed(x))## $mean.trimmed
## [1] 50.14574
##
## $sd.trimmed
## [1] 1.508397print(paste("Moyenne classique",mean(x)))## [1] "Moyenne classique 55.521587552483"print(paste("Moyenne tronquée",stat.trimmed(x)$mean.trimmed))## [1] "Moyenne tronquée 50.1457440590004"print(paste("Ecart type classique",sd(x)))## [1] "Ecart type classique 54.7485544087614"print(paste("Ecart type tronquée",stat.trimmed(x,trim=0.02)$sd.trimmed))## [1] "Ecart type tronquée 1.64596360926422"Calcul par groupe
Simulons deux groupes, un de 45 femmes et l’autre de 60 hommes, et des tailles associées à chaque individu
- La taille des femmes suit \(\mathcal N(\mu=168,\sigma=5)\)
- La taille des hommes suit \(\mathcal N(\mu=178,\sigma=6)\)
sex<-rep(c("H","F"),c(60,45))
# Forçons le type
sex<-as.factor(sex)
Tailles<-c(rnorm(60,mean=178,sd=6),rnorm(45,168,5))
# la fonction tapply applique une fonction (ici mean) à un vecteur (Tailles ) suivant les groupes définis par un facteur (ici sex)
tapply(Tailles,sex,mean)## F H
## 168.2313 178.1046tapply(Tailles,sex,function(group){max(group) - min(group)})## F H
## 16.48933 23.60910Matrices
Soit un ensemble de bien immobiliers décrits par la surface habitable et la taille du jardin.
Le prix moyen de vente d’un bien est une combinaison linéaire de la surface habitable et de la taille du jardin.
Ecrire une fonction qui prenne une matrice décrivant 100 bien immobiliers et calculer leur prix moyen par une multiplication matricielle.
Surface.maison<-rnorm(100,110,20)
Jardin<-c(rep(0,10),rnorm(90,500,100))
print(Bien.immobilier<-matrix(c(Surface.maison,Jardin),100,2))## [,1] [,2]
## [1,] 109.55591 0.0000
## [2,] 110.57813 0.0000
## [3,] 94.01510 0.0000
## [4,] 138.63233 0.0000
## [5,] 138.14143 0.0000
## [6,] 89.48352 0.0000
## [7,] 136.78644 0.0000
## [8,] 93.85385 0.0000
## [9,] 118.16659 0.0000
## [10,] 113.65089 0.0000
## [11,] 109.06524 513.5764
## [12,] 92.21736 504.4645
## [13,] 99.44175 413.0441
## [14,] 69.33212 460.9178
## [15,] 127.31404 368.4948
## [16,] 115.00425 594.7622
## [17,] 104.97474 668.0746
## [18,] 135.08181 618.7175
## [19,] 117.04207 666.6394
## [20,] 79.18027 583.2785
## [21,] 106.93312 615.5159
## [22,] 100.73174 534.9822
## [23,] 121.42528 476.9347
## [24,] 90.85923 666.3922
## [25,] 106.34573 527.8002
## [26,] 118.24440 432.8528
## [27,] 112.88979 607.7100
## [28,] 147.89230 571.5411
## [29,] 49.78234 395.5226
## [30,] 141.79390 583.3987
## [31,] 138.43799 636.5777
## [32,] 121.69063 562.2194
## [33,] 119.56118 464.7763
## [34,] 123.16169 497.4678
## [35,] 121.46915 393.0092
## [36,] 92.47942 515.8134
## [37,] 106.24149 595.6961
## [38,] 84.00115 629.2791
## [39,] 134.86774 470.8625
## [40,] 119.19671 710.7394
## [41,] 112.25967 593.3191
## [42,] 106.97329 307.8354
## [43,] 122.82767 534.0987
## [44,] 107.60899 613.2202
## [45,] 115.40730 518.3522
## [46,] 118.64515 308.2042
## [47,] 131.54673 703.0274
## [48,] 121.01429 552.5823
## [49,] 123.51015 485.0440
## [50,] 111.30326 579.5047
## [51,] 150.31160 490.3473
## [52,] 97.84190 533.1480
## [53,] 101.12303 647.8214
## [54,] 116.24790 438.6747
## [55,] 106.34961 461.7529
## [56,] 138.60534 559.8673
## [57,] 125.66135 432.3784
## [58,] 120.13451 543.3732
## [59,] 95.45837 581.0207
## [60,] 115.33757 377.5706
## [61,] 149.25973 471.0123
## [62,] 69.74540 721.5351
## [63,] 73.67174 511.2828
## [64,] 147.74275 535.5185
## [65,] 103.60451 536.9294
## [66,] 125.17250 420.0502
## [67,] 101.73469 301.2306
## [68,] 114.10435 519.6208
## [69,] 145.24085 533.4517
## [70,] 75.97805 563.8308
## [71,] 97.49121 428.5237
## [72,] 117.63067 523.7128
## [73,] 106.73060 595.4779
## [74,] 116.66114 585.5748
## [75,] 99.65090 462.3068
## [76,] 128.60115 581.5218
## [77,] 109.68172 528.2040
## [78,] 109.98603 509.9327
## [79,] 136.23954 459.7015
## [80,] 104.81777 546.6626
## [81,] 74.96266 394.1976
## [82,] 114.63110 572.6491
## [83,] 131.19909 599.6292
## [84,] 116.64220 473.5984
## [85,] 142.76511 388.6391
## [86,] 124.41127 603.5130
## [87,] 132.94598 595.9933
## [88,] 108.35875 570.0943
## [89,] 95.27601 448.7260
## [90,] 108.19926 524.9240
## [91,] 126.62309 415.8677
## [92,] 134.96509 500.9503
## [93,] 100.30142 478.5400
## [94,] 131.49495 374.2904
## [95,] 106.63464 515.7713
## [96,] 83.21664 500.2164
## [97,] 124.72148 729.0838
## [98,] 117.14646 406.1160
## [99,] 124.90158 559.8010
## [100,] 102.01682 549.4544colnames(Bien.immobilier)<-c("Surface.habitable","Taille.jardin")
Formule.regression<-c(4000,100)
print(Prix<-Bien.immobilier%*% cbind(Formule.regression))## Formule.regression
## [1,] 438223.6
## [2,] 442312.5
## [3,] 376060.4
## [4,] 554529.3
## [5,] 552565.7
## [6,] 357934.1
## [7,] 547145.8
## [8,] 375415.4
## [9,] 472666.4
## [10,] 454603.6
## [11,] 487618.6
## [12,] 419315.9
## [13,] 439071.4
## [14,] 323420.3
## [15,] 546105.7
## [16,] 519493.2
## [17,] 486706.4
## [18,] 602199.0
## [19,] 534832.2
## [20,] 375049.0
## [21,] 489284.1
## [22,] 456425.2
## [23,] 533394.6
## [24,] 430076.2
## [25,] 478162.9
## [26,] 516262.9
## [27,] 512330.2
## [28,] 648723.3
## [29,] 238681.6
## [30,] 625515.5
## [31,] 617409.7
## [32,] 542984.5
## [33,] 524722.4
## [34,] 542393.5
## [35,] 525177.5
## [36,] 421499.0
## [37,] 484535.6
## [38,] 398932.5
## [39,] 586557.2
## [40,] 547860.8
## [41,] 508370.6
## [42,] 458676.7
## [43,] 544720.5
## [44,] 491758.0
## [45,] 513464.4
## [46,] 505401.0
## [47,] 596489.7
## [48,] 539315.4
## [49,] 542545.0
## [50,] 503163.5
## [51,] 650281.1
## [52,] 444682.4
## [53,] 469274.3
## [54,] 508859.1
## [55,] 471573.7
## [56,] 610408.1
## [57,] 545883.2
## [58,] 534875.4
## [59,] 439935.5
## [60,] 499107.3
## [61,] 644140.1
## [62,] 351135.1
## [63,] 345815.2
## [64,] 644522.9
## [65,] 468111.0
## [66,] 542695.0
## [67,] 437061.8
## [68,] 508379.5
## [69,] 634308.6
## [70,] 360295.3
## [71,] 432817.2
## [72,] 522894.0
## [73,] 486470.2
## [74,] 525202.1
## [75,] 444834.3
## [76,] 572556.8
## [77,] 491547.3
## [78,] 490937.4
## [79,] 590928.3
## [80,] 473937.3
## [81,] 339270.4
## [82,] 515789.3
## [83,] 584759.3
## [84,] 513928.6
## [85,] 609924.4
## [86,] 557996.4
## [87,] 591383.3
## [88,] 490444.4
## [89,] 425976.6
## [90,] 485289.4
## [91,] 548079.1
## [92,] 589955.4
## [93,] 449059.7
## [94,] 563408.8
## [95,] 478115.7
## [96,] 382888.2
## [97,] 571794.3
## [98,] 509197.4
## [99,] 555586.4
## [100,] 463012.7hist(Prix)
Estimateur<-function(X,y){
solve(t(X)%*%X)%*%t(X)%*%y
}
print(Estimateur(X=Bien.immobilier,y=Prix))## Formule.regression
## Surface.habitable 4000
## Taille.jardin 100Listes
Accès aux éléments d’une liste
maliste<-list(A=matrix(0,2,2),B="Christophe",C=TRUE)
str(maliste)## List of 3
## $ A: num [1:2, 1:2] 0 0 0 0
## $ B: chr "Christophe"
## $ C: logi TRUEstr(maliste[1])## List of 1
## $ A: num [1:2, 1:2] 0 0 0 0str(maliste[[1]])## num [1:2, 1:2] 0 0 0 0Créer une liste d’échantillons de tailles de nourrissons:
- Chaque case de la liste est un échantillon de taille différente.
- La taille de l’échantillon suit une loi de poisson de paramètres 150.
- La liste contient 22 éléments.
# Simulons les tailles d'échantillons
taille.echantillon<-rpois(22,lambda=150)
Hopitaux<-list()
for (i in 1:22) Hopitaux[[i]]<-rnorm(taille.echantillon[i],50,2)
lapply(Hopitaux, min)## [[1]]
## [1] 45.13135
##
## [[2]]
## [1] 45.37916
##
## [[3]]
## [1] 44.04215
##
## [[4]]
## [1] 43.34796
##
## [[5]]
## [1] 44.74898
##
## [[6]]
## [1] 43.9133
##
## [[7]]
## [1] 44.03992
##
## [[8]]
## [1] 44.6237
##
## [[9]]
## [1] 44.751
##
## [[10]]
## [1] 42.24569
##
## [[11]]
## [1] 44.84959
##
## [[12]]
## [1] 44.76095
##
## [[13]]
## [1] 44.48031
##
## [[14]]
## [1] 43.7116
##
## [[15]]
## [1] 44.10751
##
## [[16]]
## [1] 45.09326
##
## [[17]]
## [1] 45.12667
##
## [[18]]
## [1] 46.08784
##
## [[19]]
## [1] 45.28996
##
## [[20]]
## [1] 45.25966
##
## [[21]]
## [1] 45.72968
##
## [[22]]
## [1] 44.33984unlist(lapply(Hopitaux, function(x) max(x)-min(x)))## [1] 10.645082 10.426090 13.596182 11.135801 10.418550 11.013927 11.462580
## [8] 10.323290 10.618242 14.139971 9.643210 10.543616 10.701138 11.057027
## [15] 11.591177 8.962536 10.982516 9.794923 10.725311 9.273790 9.222000
## [22] 10.407644for (i in 1:length(Hopitaux)) min(Hopitaux[[i]])Fonction avec plusieurs sorties
stat.des<-function(x){
return(list(moyenne=mean(x),ecart.type=sd(x)))
}
Tailles<-rnorm(70,50,2)
resultat<-stat.des(Tailles)
stat.trimmed<-function(x,trim=0.1){
n<-length(x)
k<-round(trim*n)
mean.trim<-mean(sort(x)[(k+1):(n-k)])
sd.trim<-sd(sort(x)[(k+1):(n-k)])
return(list(mean=mean.trim,sd=sd.trim))
}
echantillon<- c(rnorm(70,50,2),500,600)
stat.trimmed(echantillon)## $mean
## [1] 50.02518
##
## $sd
## [1] 1.158833echantillon["case 25"]## [1] NAVariable qualitative ordinale
Calcul des fréquences relatives et cumulées sur la vitesse dans le
jeu de données cars
data(cars)
X<-cars$speed
# Les modalités: epsilonk
epsilonk<-unique(X)
# Les effectifs: nk
nk<-table(X)
plot(nk)
# Les fréquences: fk
fk<-nk/sum(nk)
# Les fréquences cumulées: Fk
Fk<-cumsum(fk)
par(mfrow=c(1,2))
plot(Fk)
plot(ecdf(X))
Avec ggplot
library(ggplot2)
ggplot(data = cars,aes(x=dist))+geom_bar()
Quelle proportion d’une variable alétoire normale est définie par l’intervalle \([Q_1-1.5 IQR, Q_3+1.5 IQR]\) ?
Q1<-qnorm(0.25)
Q3<-qnorm(0.75)
IQR<-Q3-Q1
m.inf<- Q1-1.5*IQR
m.sup<- Q3+1.5*IQR
print(pnorm(m.sup)-pnorm(m.inf))## [1] 0.9930234# Boxplot pour distribution normale (symétrique)
X<-rnorm(1000)
boxplot(X)
# Boxplot pour distribution log normale
X<-rlnorm(1000)
ggplot(data=data.frame(salaire=X),aes(x=salaire))+
geom_boxplot(outlier.colour = "red")
Exemple du jeu de données des iris
data(iris)Graphes quantiles quantiles
x<-rnorm(50,mean=2,sd=1)
y<-rnorm(5)
qqplot(x,y)
abline(-2,1)
Estimation Covariance entre deux échantillons appariés
X<-iris[]
X<-iris[,1:4]
# Estimation par R
cov(X)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 0.6856935 -0.0424340 1.2743154 0.5162707
## Sepal.Width -0.0424340 0.1899794 -0.3296564 -0.1216394
## Petal.Length 1.2743154 -0.3296564 3.1162779 1.2956094
## Petal.Width 0.5162707 -0.1216394 1.2956094 0.5810063X<-scale(X,center=TRUE,scale=FALSE)
n<-nrow(X)
# Notre estimation
1/(n-1)*t(X)%*%X## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 0.6856935 -0.0424340 1.2743154 0.5162707
## Sepal.Width -0.0424340 0.1899794 -0.3296564 -0.1216394
## Petal.Length 1.2743154 -0.3296564 3.1162779 1.2956094
## Petal.Width 0.5162707 -0.1216394 1.2956094 0.5810063cor(X)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 1.0000000 -0.1175698 0.8717538 0.8179411
## Sepal.Width -0.1175698 1.0000000 -0.4284401 -0.3661259
## Petal.Length 0.8717538 -0.4284401 1.0000000 0.9628654
## Petal.Width 0.8179411 -0.3661259 0.9628654 1.0000000Graphe de dispersion (biplot)
library(MASS)
data(crabs)
help(crabs)
pairs(crabs)
pairs(crabs[,4:8])
cor(crabs[,4:8])->correl
colMeans(correl)## FL RW CL CW BD
## 0.9676825 0.9178676 0.9699622 0.9656384 0.9655696(colSums(correl)-1)/4## FL RW CL CW BD
## 0.9596031 0.8973345 0.9624528 0.9570480 0.9569620Faisons un changement de variables pour raisonner en relatif par
rapport à la variable longueur de carapace. Cela permet de
se focaliser sur les proportions des mesures corporelles.
crabs.new<-(crabs[,4:8]/crabs$CL)[,-3]
names(crabs.new)<-paste(names(crabs.new),"/CL",sep="")
group<-as.factor(paste(crabs$sex,crabs$sp,sep="-"))
pairs(crabs.new,col=c("cyan","orange","blue","orange3")[group])