Simulation des données

  1. Commencer par simuler \(n=500\) variables aléatoires suivant un mouvement brownien standard calculé en \(p=300\) points de discrétisation. On les notera \(X_1,...,X_n\).
MvB=function(n,p)
{
  # Retourne une matrice de taille n*p (X_i(t_j)) i=1,...,n, j=1,...,p
  # tel que X_i est un mouvement brownien standard.
  X=rnorm(n*p)/sqrt(p-1)
  X=matrix(X,ncol=p)
  MvB=matrix(0,ncol=p,nrow=n)
  for (i in 1:n)
  {
    MvB[i,]=cumsum(X[i,])
  }
  MvB
}
n = 500
p = 300
X = MvB(n,p)
t = seq(0,1,length.out = p)
plot(c(0,1),range(X),type='n',xlab='t',ylab=expression(X[i](t)),main='n mouvements browniens sur [0,1]')
for (i in 1:n){
  points(t,X[i,],type='l',col=i)
}

  1. Pour \(\beta^*(t) = \sin(10t)\), générer des variables aléatoires \(Y_1,...,Y_n\) à partir de \(X_1,...,X_n\) telles que \[ Y_i = \int_0^1X_i(t)\beta^*(t)dt+\varepsilon_i, \] où \(\varepsilon_1,...,\varepsilon_n\sim_{i.i.d}\mathcal N(0,\sigma^2)\) avec \(\sigma=0.1\) et \(\varepsilon_1,...,\varepsilon_n\) indépendantes de \(X_1,..., X_n\). On pourra, lorsque \(t_1,...t_p\) est une grille régulière de \([0,1]\) formée de \(p\) points, approcher l’intégrale par une somme \[ \int_0^1 f(t)dt \approx \frac1p\sum_{j=1}^pf(t_j). \]
sigma = 0.1
epsilon = rnorm(n,0,sigma)
betas = function(t){sin(10*t)}
Y  = X%*%betas(t)/p + epsilon
  1. Calculer et tracer l’estimateur \(\widehat\beta_m\) pour différentes valeurs de \(m\) (par exemple \(m=1\), \(m=5\), \(m=31\)). \[ \widehat\beta_m = {\arg\min}_{\beta\in S_m}\gamma_n(\beta) \text{ où }\gamma_n(\beta)=\left\{\frac1n\sum_{i=1}^n(Y_i-\langle \beta,X_i\rangle)^2\right\} \] où \(S_m={\rm Vect}\{\varphi_1,...,\varphi_m\}\) est la base de Fourier définie au TP2. On pourra reprendre les codes du TP2 en les adaptant.
mmax = 100
Dmax = 2*mmax+1
phiX = matrix(NA,n,Dmax) # matrice contenant <X_i,phi_k>
phi = matrix(1,Dmax,p) # matrice contenant phi_k(t_l)
phiX[,1] = rowMeans(X)
for (j in 1:mmax){
  phiX[,2*j] = (sqrt(2)/p)*X%*%cos(2*pi*j*t)
  phiX[,2*j+1] = (sqrt(2)/p)*X%*%sin(2*pi*j*t)
  phi[2*j,] = sqrt(2)*cos(2*pi*j*t)
  phi[2*j+1,] = sqrt(2)*sin(2*pi*j*t)
}
G = crossprod(phiX)/n
hatb = crossprod(phiX,Y)/n
m1 = 1
hata1 = hatb[1]/G[1,1]
hatbeta1 = rep(hata1,p)
m2 = 5
hata5 = solve(G[1:m2,1:m2],hatb[1:m2])
hatbeta5 = hata5%*%phi[1:m2,]
m3 = 31
hata31 = solve(G[1:m3,1:m3],hatb[1:m3])
hatbeta31 = hata31%*%phi[1:m3,]
plot(t,betas(t),col="darkgreen",type='l',lwd=2,lty=2)
abline(h=hata1, col='blue')
points(t,hatbeta5,type='l',col="green")
points(t,hatbeta31,type='l',col="red")
legend('bottomleft',c(paste(expression(beta),'*'),expression(hat(beta)[1]),expression(hat(beta)[5]),expression(hat(beta)[31])),lty=c(2,rep(1,3)),lwd=c(2,rep(1,3)),col=c('darkgreen','blue','green','red'))

  1. Tracer, pour \(\kappa = 2.5\), les fonctions \[ m\mapsto \gamma_n(\widehat\beta_m)+\kappa\frac{\sigma^2m}n \] et

\[ m\mapsto \gamma_n(\widehat\beta_m)\left(1+\kappa\frac{m}n\right). \] Puis tracer les estimateurs obtenus en sélectionnant \(m\) par les deux critères de sélection de modèle (variance connue et inconnue).

kappa = 2.5
gamma = rep(0,Dmax)
for (m in 1:Dmax){
  hatam = solve(G[1:m,1:m],hatb[1:m])
  hatbetam = hatam%*%phi[1:m,]
  gamma[m] = mean((Y-tcrossprod(X,hatbetam)/p)^2)
}
plot(gamma,col='blue',type='b',pch=16,xlab='m',ylab='crit(m)')
points(gamma+kappa*sigma^2*(1:Dmax)/n,col='green',type='l')
points(gamma*(1+kappa*(1:Dmax)/n),col='orange',type='l')
legend('topright',c(expression(gamma[n](hat(beta)[m])),'crit1','crit2'),lty=rep(1,3),col=c('blue','green','orange'))

hatm1 = which.min(gamma+kappa*sigma^2*(1:Dmax)/n)
hatm1
## [1] 6
hatm2 = which.min(gamma*(1+kappa*(1:Dmax)/n))
hatm2
## [1] 6
hatahatm1 = solve(G[1:hatm1,1:hatm1],hatb[1:hatm1])
hatbetahatm1 = hatahatm1%*%phi[1:hatm1,]
hatahatm2 = solve(G[1:hatm2,1:hatm2],hatb[1:hatm2])
hatbetahatm2 = hatahatm2%*%phi[1:hatm2,]
plot(t,betas(t),col="darkgreen",type='l',lwd=2,lty=2)
points(t,hatbetahatm1,type='l',col="green",lwd=2)
points(t,hatbetahatm2,type='l',col="orange")
legend('bottomleft',c(expression(beta),expression(hat(beta)[hat(m)[1]]),expression(hat(beta)[hat(m)[2]])),lty=c(2,1,1),lwd=c(2,2,1),col=c('darkgreen','green','orange'))

Application à l’évaluation de la concentration en nitrate des eaux usées

Nous allons étudier appliquer la méthode d’estimation à l’a préd’estimation de la concentration en nitrates des eaux usées. Ces données sont celles étudiées dans l’article

H.N. Pham et al. Estimation simultanée et en ligne de nitrates et nitrites par identification spectrale UV en traitement des eaux usées, L’eau, l’industrie, les nuisances, n. 335, p. 61-69, 2010.

Nous commençons par charger et centrer les données. Les données concernant la variable à prédire sont enregistrées dans la première colonne du fichier Y.dat.

Xnc <- read.table("X.dat")
Ync <- read.table("Y.dat")
Ync <- Ync[,1]
n <- length(Ync)

# Representation des donnees selon la valeur de Y
palette(rainbow(14,start=0,end=4/6))
plot(c(0,1),range(Xnc),type='n',xlab='wavelength',ylab='absorbance')
for (i in 1:n){
  points(seq(0,1,length.out=ncol(Xnc)),Xnc[i,],type='l',col=floor(Ync[i]))
}
ry = range(Ync)
legend("topright",paste(floor(ry[1]):floor(ry[2]),'<=Y<',floor(ry[1]+1):floor(ry[2]+1)),col=1:9,lty=rep(1,9))

## Centrage des données
Xbar <- colMeans(Xnc)
points(1:ncol(Xnc),Xbar,col=2,type='l')

X <- scale(Xnc,scale=F)
plot(c(1,ncol(Xnc)),range(X),type='n',xlab='wavelength',ylab='',main="Données centrées")
for (i in 1:n){
  points(1:ncol(Xnc),X[i,],type='l',col=floor(Ync[i]))
}

Y1 <- Ync-mean(Ync)

Appliquer la méthode d’estimation de l’exercice précédent et proposer une méthode de prédiction de la variable Y.