# *********************************************************************
# **  Licences Sciences pour la Santé - L2 - Statistiques - TD1      **
# **             Simulation de lois Student et Fisher		         **
# *********************************************************************


rm(list=ls())

# chemin du répertoire
Dir <- "U:/enseignement/SciencesSante/L2/TD/TD1_Simu_loi_Student_Fisher/"


sink(paste0(Dir,"TD1.lis"))

date()
system("hostname", intern=T)

print(sessionInfo())


#----------------------------------------------------------------------------------------------------------------------------------------
## 1 ==>  Simulation d'une loi normale centrée réduite

# simulation d'1 valeur + histogramme
x = rnorm(1,mean=0,sd=1)
hist(x)

# simulation de 10 valeurs + histogramme
x = rnorm(10,mean=0,sd=1)
hist(x)

# simulation de 100 valeurs + histogramme
x = rnorm(100,mean=0,sd=1)
hist(x)

# simulation de 1000 valeurs + histogramme
x = rnorm(1000,mean=0,sd=1)
hist(x)

# simulation de 1 000 000 valeurs + histogramme
x = rnorm(10^6,mean=0,sd=1)
hist(x)


#------- comparer la simulation avec la densité théorique

# moyenne

mean(x) # espérance théorique = 0

# variance

var(x) # variance théorique = 1

# distribution

hist(x,freq=FALSE)

nx <- seq(-100,100,length=5000)
lines(nx,dnorm(nx,mean=0,sd=1),col="red")

# quantiles

qqnorm(x, pch = 1, frame = FALSE)
qqline(x, col = "steelblue", lwd = 2)


#---- faire son propre graph quantile-quantile

x = rnorm(10,mean=0,sd=1)

x = sort(x)

x

# nombre de valeurs

nb <- length(x)

nb

# rang des valeurs

rang <- 1:nb

rang

# percentiles

perc <- (rang-0.5)/nb

perc

# Z-scores : valeurs d'une loi normale centrée réduite correspondant aux percentiles calculés

Zscore <- qnorm(perc,mean=0,sd=1)


# QQplot

par(mfrow=c(1,2))
plot(Zscore,x,xlab="théorique",ylab="observé")


qqnorm(x, pch = 1, frame = FALSE)
qqline(x, col = "steelblue", lwd = 2)



#----------------------------------------------------------------------------------------------------------------------------------------
## 2. ==>  Simulation d'une loi du khi deux

# Loi du khi-deux à 1 degré de liberté

x2 <- x^2


mean(x2) # espérance théorique = df

var(x2) # variance théorique = 2*df


hist(x2,freq=FALSE)

nx <- seq(-100,100,length=5000)
lines(nx,dchisq(nx,df=1),col="red")


# Loi du khi-deux à 3 degrés de liberté

n = 1000

x1 <- rnorm(n,mean=0,sd=1)
x2 <- rnorm(n,mean=0,sd=1)
x3 <- rnorm(n,mean=0,sd=1)

x <- x1^2 + x2^2 + x3^2



mean(x) # espérance théorique = df

var(x) # variance théorique = 2*df


hist(x,freq=FALSE)

nx <- seq(-100,100,length=5000)
lines(nx,dchisq(nx,df=3),col="red")



# qqplot
par(mfrow=c(1,1))
y <- qchisq(ppoints(length(x)),df=3)
qqplot(x,y)


#----------------------------------------------------------------------------------------------------------------------------------------
## 3. ==>  Simulation d'une loi de Student

n=10^5

df1 = 50



x <- rnorm(n,mean=0,sd=1)/sqrt(rchisq(n,df=df1)/df1)



# moyenne

mean(x) # espérance théorique 0


# variance

var(x) # variance théorique (df/(df-2))

df1/(df1-2)



# distribution

hist(x,freq=FALSE)

nx <- seq(-15,15,length=5000)
lines(nx,dt(nx,df=3),col="red")

lines(nx,dnorm(nx,mean=0,sd=df1/(df1-2)),col="blue3")


# qqplot
par(mfrow=c(1,1))
y <- qt(ppoints(length(x)),df=df1)
qqplot(x,y)




#----------------------------------------------------------------------------------------------------------------------------------------
## 4. ==>  Simulation d'une loi de Fisher


n=10^5

df1 = 10
df2 = 15


x <- rchisq(n,df=df1)/df1/(rchisq(n,df=df2)/df2)



# moyenne

mean(x) # espérance théorique (df2/(df2-2))

df2/(df2-2)



# variance

var(x) # variance théorique (2*df2^2*(df1+df2-2)/(df1*(df2-2)^2*(df2-4)))

2*df2^2*(df1+df2-2)/(df1*(df2-2)^2*(df2-4))




# distribution

hist(x,freq=FALSE)

nx <- seq(0,30,length=5000)
lines(nx,df(nx,df1=df1,df2=df2),col="red")




# qqplot
par(mfrow=c(1,1))
y <- qf(ppoints(length(x)),df1=df1,df2=df2)
qqplot(x,y)




date()

sink()