TP3

On change le chemin courrant pour se placer là où on a sauvegarder le fichier polynomesorthogonaux

setwd("E:/TP Monte Carlo")
source("TP3_codepolynomesorthogonaux.R")

Exercice 1.3

\(\mathbb E(Y|X) = f(X) = \dfrac{\max(1,X)-0.5}{2}\)

N=20000 # nombre de simulations

X=rexp(N,1)
Y=numeric()
for (i in seq(1,N)){    
    Y[i]=(max(X[i],1)-1/2)*runif(1)
}

Phi=cbind(rep(1,N), R1(X), R2(X),R3(X), R4(X), R5(X))
colnames(Phi)=rep("",6)

A=t(Phi)%*%Phi
B=t(Phi)%*%Y
alpha=solve(A,B)
alpha

##         [,1]
##   0.43441597
##  -0.38114122
##   0.12031715
##   0.02889540
##   0.04270572
##   0.05170193

Espcond=Phi%*%alpha
hat.f<-function(x){alpha[1]+R1(x)*alpha[2]+R2(x)*alpha[3]+R3(x)*alpha[4]+R4(x)*alpha[5]+R4(x)*alpha[5]}

x=sort(X)
y=1/2*(pmax(1,x)-1/2)
plot(X,Espcond,col="blue",type="p",pch=".")
lines(x,y,col="red",type="l")

approximation d’une fonction non dérivable par des polynômes…

plot(X,Y)
plot(hat.f,xlim=c(0,8),add=TRUE,col="blue")

#lines(x,y,type="l",col="red")

Exercice 1.4

\(U\sim \gamma (2,1)\) on a \[F(u) = 1 - (1+u)\exp(-u)\] Ainsi \(F^{-1}(\mathcal U)\) solution en \(u\) de \((1+u)\exp(-u) - 1 - \mathcal U = 0\), où \(\mathcal U\) uniforme sur \([0,1]\) et \(F^{-1}(\mathcal U)\) à même loi que \(U\).

N=10000 # nombre de simulations
X=runif(N);
U=rep(0,N);

f<-function(t,x){1-(x+1)*exp(-x)-t}
fx<-function(t){uniroot(function(x){f(t,x)},c(0,20))$root}
U=sapply(runif(N),fx)

Vérifions que \(U\) suit la loi gamma \(\gamma(2,1)\)

hist(U,freq=FALSE)
plot(function(x){dgamma(x,2,1)},col="red",add=TRUE,xlim=c(0,10))

Y=X+U

# Base classique

Phi=cbind(rep(1,N), X, X^2, X^3,X^4)
colnames(Phi)=rep("",5)

A=t(Phi)%*%Phi
B=t(Phi)%*%Y
alpha1=solve(A,B)
alpha1

##         [,1]
##   2.06584749
##   0.48126778
##   0.76802639
##  -0.28197383
##  -0.00807619

hat.f1=function(x){alpha1[1]+x*alpha1[2]+x^2*alpha1[3]+x^3*alpha1[4]+x^4*alpha1[5]}


Espcond1=Phi%*%alpha1;
mean(Espcond1)

## [1] 2.48894

###On verifie que E(E(Y|X)) = 2+E(X) = 2.5

x=sort(X)
plot(X,Espcond1,type="p")
lines(x,x+2,col="red")

Phi=cbind(rep(1,N), R1(X), R2(X), R3(X),R4(X))
colnames(Phi)=rep("",5)

A=t(Phi)%*%Phi
B=t(Phi)%*%Y
alpha2=solve(A,B)
alpha2

##        [,1]
##   2.1974970
##   2.2974679
##  -4.7024442
##   2.4671547
##  -0.1938279

hat.f2=function(x){alpha2[1]+R1(x)*alpha2[2]+R2(x)*alpha2[3]+R3(x)*alpha2[4]+R4(x)*alpha2[5]}


Espcond2=Phi%*%alpha2;
mean(Espcond2)

## [1] 2.48894

x=sort(X)
plot(X,Espcond2,type="p")
lines(x,x+2,col="red")

plot(X,Y,pch=".")
plot(hat.f1,add=TRUE,col="red")
plot(hat.f2,add=TRUE,col="blue")