######################################################################
#### Estimation de densité
###############################################################

## Partie A : simulation des données
n=10^3
x <- seq(0,5, l=2000)

## question 1
Y <- rnorm(n, mean=2.5)
dY <- dnorm(x, mean=2.5, sd=1)
plot(x,dY,t='l', col='red')

## question 2-a)
Mu <- c(1,3,4)
Sd <- c(0.5,0.3,0.2)
dX <- 0.2*dnorm(x, mean=Mu[1], sd=Sd[1]) +0.5*dnorm(x, mean=Mu[2], sd=Sd[2])+0.3*dnorm(x,mean=Mu[3], sd=Sd[3])
plot(x,dX,t='l', col='red')

## question 2-b)

U <- sample(x=c(1,2,3), size=n, replace=TRUE, prob=c(0.2,0.5,0.3))
X <- vector(l=n)
for(i in 1:n) X[i] <- rnorm(1, mean=Mu[U[i]], sd=Sd[U[i]])


## Partie B : estimation par histogrammes

## question 1-a)
par(mfrow=c(2,2))
for(D in c(5,15, 50, 300)){
		hist(Y,D, freq=F)
	lines(x,dY,t='l', col='red')
	}
	
## question 2)
par(mfrow=c(2,2))
for(D in c(5,15, 50, 300)){
		hist(Y,D, freq=F)
	lines(x,dY,t='l', col='red')
	}
	
## question 4)
library(histogram)
par(mfrow=c(1,2))
histogram(X, type='irregular', verbose=F)
lines(x,dX,t='l', col='red')
				
histogram(Y, type='irregular', verbose=F)
lines(x,dY,t='l', col='red')

## Partie C : estimation par noyaux

## question 1-a)
kernX <- density(X)
plot(kernX)
lines(x,dX,t='l', col='red')

## question 1-b)
kernX$bw

## question 1-c)
bw.nrd0(X)
bw.nrd(X)
bw.bcv(X) 
bw.ucv(X)
bw.SJ(X)

## question 1-d)
par(mfrow=c(2,2))
for(b in c(0.01,0.1, 0.22, 2)){
		plot(density(X,bw=b))
	lines(x,dX,t='l', col='red')
	}

## question 2)
kernY <- density(Y)
plot(kernY)
lines(x,dY,t='l', col='red')
kernY$bw
bw.nrd0(Y)
bw.nrd(Y)
bw.bcv(Y) 
bw.ucv(Y)
bw.SJ(Y)

par(mfrow=c(2,2))
for(b in c(0.02,0.08, 0.2, 2)){
		plot(density(Y,bw=b))
	lines(x,dY,t='l', col='red')
	}

## question 3)
N=10^5
P <- sample(x=c(1,2,3), size=N, replace=TRUE, prob=c(0.2,0.5,0.3))
X1 <- vector(l=N)
#Mu <- c(mu1,mu2,mu3)
#Sd <- c(sd1,sd2,sd3)
for(i in 1:N) X1[i] <- rnorm(1, mean=Mu[P[i]], sd=Sd[P[i]])
	
par(mfrow=c(2,2))
for(h in c(0.05,0.2,0.4,0.8)){
	plot(x,dX, col='red', t='l')
	d <- density(X1, bw=h)
	lines(d)
	#lines(x, dX, col='red')
	
}

		
## Partie D

## question 1
library(datasets)	
?precip
plot(density(precip))
hist(precip,15)	

## question 2
library(MASS)
?Aids2
plot(density(Aids2$age))
hist(Aids2$age,20)


######################################################################
#### Tests NP : cas pratiques 
###############################################################



library(MASS)
library(nortest)

#########  Exercice 1 ################

## poids avant traitement
X <- anorexia[,2]
## poids après traitement
Y <- anorexia[,3]

# question 1-a
qqnorm(X-Y)
qqline(X-Y)
lillie.test(X-Y)
t.test(X-Y)

## question 1-b
wilcox.test(X-Y)

##  Difference de poids avant après traitement en fonction du traitement
Z <- anorexia[,3] -anorexia[,2]
Z1 <- Z[(anorexia$Treat=='Cont')]
# alternative : Z1 <- Z[which(anorexia$Treat=='Cont')]
Z2 <- Z[(anorexia$Treat=='FT')]

## Question 2-a
qqplot(Z1,Z2)
abline(0,1)

## Question 2-b

lillie.test(Z1)
lillie.test(Z2)
t.test(Z1,Z2)

## Question 2-c 
wilcox.test(Z1,Z2)

### Combien y a-t-il d'ex aequos?
length(unique(Z1))
length(Z1)
length(unique(Z2))
length(Z2)


#########  Exercice 2 ################

## Différence entre taux de criminalité chez les 18-24, entre états du sud et du nord
head(UScrime)

X <- UScrime$Prob[(UScrime$So==1)]
Y <- UScrime$Prob[(UScrime$So==0)]

## Intuition a partir de qqplot 
qqplot(X,Y)
abline(0,1)
## D'après qqplot(X,Y) : différence de médiane significative 

## test de Mann Whitney

wilcox.test(X,Y)


#########  Exercice 3 ################


head(Melanoma)
X <- Melanoma$time[Melanoma$sex==0]
Y <- Melanoma$time[Melanoma$sex==1]

## Intuition
qqplot(X,Y)
abline(0,1)
##  test de mann whitney 
wilcox.test(X,Y)


#########  Exercice 4 ################

head(Animals)
X <- Animals$body
Y <- Animals$brain
plot(X,Y)

## Question a 
mod <- lm(Y~X)
plot(X,residuals(mod))
qqnorm(residuals(mod))
qqline(residuals(mod))

## Question b 
cor.test(X,Y, method='spearman', alternative='greater')
length(unique(X))
length(X)
length(unique(Y))
length(Y)