#Paramètre à ne pas changer (nécessaire à la simulation) Delta=0.01; b<-function(x) -x; sigma<-function(x) 1; T<-100; t<-seq(0,T,Delta); x=seq(0,5,0.02) X<-simulation(T) # on simule un processu d'Ornstein-Uhlenbeck sur [0,T] tt<-1:T; y<-tt;z<-y; w<-y; ww<-y; for (i in 1:T) { y[i]<-bayes(X[1:(i/Delta)],2) z[i]<-bayes(X[1:(i/Delta)],0) w[i]<-bayes(X[1:(i/Delta)],1) ww[i]<-emv(X[1:(i/Delta)]) } plot(tt,y,type='l',ylab='bayes',xlab='t',ylim=c(0,4)) lines(tt,z,col='blue',lty=2) lines(tt,w,col='darkgreen',lty=3) lines(tt,ww,col='pink',lty=5) lines(tt,2+0*tt,col='red') lines(tt,2+1.65*sqrt(8/tt),col='green') lines(tt,2-1.65*sqrt(8/tt),col='green') ## boxplot y=seq(1,100,1); z=y; w=y; ww=y; for (i in 1:100) { X=simulation(100) y[i]=bayes(X,2) z[i]=bayes(X,1) w[i]=bayes(X,0) ww[i]=emv(X) } boxplot(ww,y,z,w,ylim=c(0,4),names=c('EMV','2','1','0')) #calcul de la vraisemblance pour X en un point theta. vrais<-function(theta,X) {X1<-X[1:(length(X)-1)]; X2<-X[2:length(X)]; DX<-X2-X1; Sigma<-sigma(X1); I1<-sum(b(X1)*theta*DX/Sigma); I2<-Delta*sum(b(X1)^2*theta^2/Sigma^2); return(I1-I2/(2)) } emv<-function(X){ y<-seq(0,5,0.02); v<-y; for (i in 1:length(y)) {v[i]<-vrais(y[i],X)} e<-y[which.max(v)]; return(e) } bayes<-function(X,l){ # calcul de l'estimateur bayésien pour un à priori exponentiel, paramètre l. y<-seq(0,5,0.02); v=y; for (i in 1:length(y)) {v[i]<-vrais(y[i],X); } b<- sum(y*v* dnorm(y,l,1))/sum(v*dnorm(y,l,1)); return(b) } #Simulation (OU avec b=2) simulation<-function(T) { N<-rnorm(length(t)-1,0,sqrt(Delta)); x[1]=1; for (i in 1:(length(t)-1)) { x[i+1]=x[i]+2*b(x[i])*Delta+sigma(x[i])*N[i]; } return(x) }