-
x <- rexp(1000) # paramètre 1 par défaut
mean(x)
sd(x)
quantile(x,c(0.25,0.5,0.75))
-
{hist(x,
xlim=c(0,6), # fixer pour plusieurs graphiques
col="blue", # la couleur
freq=F, # fréquences plutôt que effectifs
main="loi exponentielle",# titre du graphique
ylab="")} # pas d'étiquette en ordonnée
d <- density(x) # estimation de densité
par(new=T) # pour superposer les graphiques
{plot(d,col="green",xlim=c(0,6),
axes=F,main="",xlab="",ylab="")}
curve(exp(-x),col="red",add=T) # ajoute la vraie densité
-
f <- ecdf(x) # fonction de répartition
plot(f,main="fonction de répartition",col="blue")
curve(1-exp(-x),col="red",add=T)
-
i <- floor(x) # parties entières
t <- table(i) # valeurs et effectifs
val <- as.numeric(row.names(t))# valeurs
fr <- as.vector(t) # effectifs
fr <- fr/sum(fr) # fréquences
pr <- dgeom(val, 1-exp(-1)) # probabilités
fp <- rbind(fr,pr) # matrice frequences-probabilites
# diagramme en barres double
barplot(fp,beside=T,col=c("blue","red"))
-
ks.test(x,"pexp",1) # appliquer le test
ks.test(x,"pexp",1)$p.value # extraire la p-valeur
-
ech <- 1000 # nombre d'échantillons
n <- 20 # taille d'échantillon
lambda <- 2 # valeur du paramètre
X <- rexp(ech*n,rate=lambda) # tirages aléatoires
X <- matrix(X,ech,n) # redimensionner
T <- 1/rowMeans(X) # inverse des moyennes
Tp <- T*(n-1)/n # estimations débiaisées
dT <- density(T) # estimation de densité
dTp <- density(Tp)
rep <- lambda*c(0.5,2) # intervalle de représentation
{plot(dT,col="blue",xlim=rep,
main="estimateurs de lambda",xlab="")}
lines(dTp,col="green",xlim=rep)
abline(v=lambda,col="red") # vraie valeur
-
ic.moy <- function(X,nc=0.95){
# retourne l'intervalle de confiance pour lambda
# basé sur la fonction t.test
#
ic <- t.test(X,conf.level=nc) # appel de la fonction
ic <- ic$conf.int # extraire l'intervalle
ic <- as.vector(ic) # conserver les deux bornes
ic <- 1/ic # inverser les deux bornes
ic <- sort(ic) # remettre dans l'ordre
return(ic) # retourner l'intervalle
}
-
ic.sym <- function(X,nc=0.95){
# retourne l'intervalle de confiance pour lambda
# basé sur la loi gamma
#
a2 <- (1-nc)/2 # probabilité des ailes
n <- length(X) # taille d'échantillon
T <- 1/mean(X) # estimation
q1 <- qgamma(a2,shape=n) # quantile inférieur
q2 <- qgamma(1-a2,shape=n) # quantile supérieur
ic <- c(q1,q2) # intervalle de quantiles
ic <- ic*T/n # intervalle de confiance
return(ic) # retourner l'intervalle
}
-
test.gamma <- function(X,la0){
# retourne deux p-valeurs pour le test
# de lambda=la0 contre lambda>la0
#
p1 <- {t.test(X, # test de la moyenne
alternative="less", # contre moyenne inférieure
mu=1/la0)$p.value} # extraire la p-valeur
T <- 1/mean(X) # estimation
S <- n*la0/T # statistique de test
p2 <- pgamma(S,shape=n) # loi sous l'hypothèse nulle
return(c(p1,p2)) # retourner les deux valeurs
}
-
ech <- 1000 # nombre d'échantillons
n <- 20 # taille d'échantillon
lambda <- 2 # valeur du paramètre
X <- rexp(ech*n,rate=lambda) # tirages aléatoires
X <- matrix(X,ech,n) # redimensionner
T <- 1/rowMeans(X) # inverse des moyennes
Tp <- T*(n-1)/n # estimations débiaisées
T1 <- (rowMeans(X^2)/2)^(-1/2) # estimation par les carrés
T2 <- log(2)/apply(X,1,median) # estimation par la médiane
A <- cbind(T,Tp,T1,T2) # faire un tableau
{boxplot(data.frame(A),names= # représenter les boxplots
c("T","Tp","T1","T2"))}
abline(h=lambda,col="red") # vraie valeur
-
lambda <- 2 # paramètre
n <- 1000 # taille des échantillons
x <- rexp(n,rate=lambda) # échantillon exponentiel
x <- sort(x) # valeurs triées
x <- x[-n] # supprimer la dernière
F <- (1:(n-1))/n # ordonnées de la fonction
y <- -log(1-F) # changement de variable
plot(x,y, pch=".", col="blue") # couples de points
reg<-lm(y~x) # régression linéaire
abline(ref,col="red") # tracer la droite
-
reg.est <- function(X){
# retourne l'estimation de lambda par
# régression non linéaire sur la fonction
# de répartition empirique de X
#
n <- length(X) # taille d'échantillon
x <- sort(X) # valeurs triées
x <- x[-n] # supprimer la dernière
F <- (1:(n-1))/n # ordonnées de la fonction
y <- -log(1-F) # changement de variable
reg <- lm(y~x) # regression linéaire
return(reg$coefficients[2]) # retourner la pente
}
compare.est <- function(lambda,E,n){
# retourne l'erreur quadratique moyenne de
# 5 estimateurs de lambda, calculée sur
# E échantillons de taille n
#
X <- rexp(E*n,rate=lambda) # tirages aléatoires
X <- matrix(X,E,n) # redimensionner
T <- 1/rowMeans(X) # inverse des moyennes
Tp <- T*(n-1)/n # estimations débiaisées
T1 <- (rowMeans(X^2)/2)^(-1/2) # estimation par les carrés
T2 <- log(2)/apply(X,1,median) # estimation par la médiane
Tr <- rep(0,E) # estimation par régression
for (i in 1:E){ # pour chaque échantillon
x <- X[i,] # extraire l'échantillon
Tr[i]<-reg.est(x) # estimation
} # fin de boucle
A <- rbind(T,Tp,T1,T2,Tr) # faire un tableau
A <- (A-lambda)^2 # carrés des erreurs
q <- sqrt(rowMeans(A)) # erreurs quadratiques
return(q) # retourner les erreurs
}
> compare.reg(2,1000,100)
T Tp T1 T2 Tr
0.2068059 0.2029628 0.2286295 0.3020214 0.2821136
-
n <- 1000 # taille des échantillons
X <- rexp(n); Y <- rexp(n) # échantillons
ind <- which(Y>((1-X)^2)/2) # indices à conserver
U <- X[ind] # valeurs conservées
hist(U^2,breaks=20,freq=F) # histogramme
{curve(dchisq(x,1),col="red", # densité de la loi du chi-deux
add=T)}
-
m <- length(U) # longueur du nouvel échantillon
S <- 2*rbinom(m,1,0.5)-1 # signes aléatoires
Z <- U*S # nouvelles valeurs
f <- ecdf(Z) # fonction de répartition
plot(f,col="blue",main="loi normale")
curve(pnorm(x),col="red",add=T)
-
qqnorm(Z)
ks.test(Z,"pnorm")
shapiro.test(Z)
-
n <- 1000 # taille des échantillons
X<-rexp(n)
Y<-rexp(n)
Z<-rexp(n)
S <- X/(X+Y+Z)
T <- (X+Y)/(X+Y+Z)
plot(S,T, pch=".", col="blue") # couples de points
ks.test(S,"pbeta",1,2)
ks.test(T,"pbeta",2,1)
-
n <- 1000; # taille des échantillons
X<-runif(2*n,0,1) # tirages
X <- matrix(X,2,n) # 2 lignes, n colonnes
X <- apply(X,2,sort) # trier chaque colonne
Sp <- X[1,] # première ligne
Tp <- X[2,] # seconde ligne
ks.test(S,Sp)
ks.test(T,Tp)
wilcox.test(S,Sp)
wilcox.test(T,Tp)
-
lambda <- 0.5 # paramètre
n <- 1000 # taille des échantillons
T <- rexp(n,rate=lambda) # échantillon exponentiel
p <- exp(-T) # paramètres
Y <- rgeom(n,p) # échantillon de géométriques
boxplot(Y) # des valeurs trop extrêmes
quantile(Y,seq(0.1,0.9,by=0.1))# déciles
max(Y) # maximum
-
tr <- 10 # valeur de troncature
Y[which(Y>tr)] <- tr # troncature
pr <- lambda*beta((1+lambda)*rep(1,tr),1:tr)
pr <- append(pr,1-sum(pr)) # probabilités théoriques
t <- table(Y) # effectifs observés
chisq.test(t,p=pr) # test du chi-deux
