############################ Fonctions annexes
# Fonction lldist(long,lat,long0,lat0) qui calcule les distances d'une ville 
# de longitude et latitude long0 et lat0 (degres) aux villes dont 

# les longitudes et latitudes sont dans les vecteurs long et lat
# Et fonction plotFr(x,y,z) qui fait un trace couleur avec les frontieres
# Valeur z[i] au point de longitude et latitude en degres (x[i],y[i]).

source("./KrigRout.R")

########################## Chargement et tracé des données
############### Table D contenant longitudes, lattitudes (degres) 
############### et temperatures (D$lon,D$lat,D$data)
############### Table Fronti contenant les frontieres.
D=read.table("./FranceTemp.txt",h=T)
Fronti=read.table("./France-borders.txt",h=T)
# Trace des frontieres
plotpoints(D$lon,D$lat,D$data)
lines(Fronti$x,Fronti$y)

################# Calcul de la matrice de distances ################# 
NT = length(D$data)
Dist=lldist(D$lon,D$lat,D$lon,D$lat)
VarT=var(D$data);
MT=mean(D$data);
T0=D$data-MT; 

################ Premier modele ##################
#Calcul de la log-vraisemblance changée de signe fois 2
cost = function(par) {
  R = ; # Matrice de variances covariances
  ll = # log-vraisemblance
  return(ll) 
}

# Estimation du modele avec optim
par0 = # Initialisation
paropt=optim(par0,cost,gr=NULL,method="L-BFGS-B",lower=??,upper=??)
parest=paropt$par;
cat("paramètre estimé =",parest)

# Vérification graphique
x=seq(parest-10,parest+10,length.out=40) 
plot(x,lapply(x,cost))

# Calcul de V0=inv(Rxx)(x-m) pour les predictions
R= # Matrice de variances covariances
V0= # 
# Calcul de la prediction m+Rxy.V0 a Rennes (longitude/latitude: -1.7 et 48):
Rxy=
cat("Température à Rennes =",) #completer

# Fabrication de la grille 100X100 pour la prédiction partout
la=seq(min(D$lat),max(D$lat),length.out=100) # les lattitudes concernees
lo=seq(min(D$lon),max(D$lon),length.out=100) # Les longitudes concernees
Glon=rep(lo,times=length(la));Glat=rep(la,each=length(lo)); # la grille en vecteur

# Caclul des predictions sur la grille et trace
tempre=Glon;
for (i in 1:length(Glon)){
  Rxy=lldist(D$lon,D$lat,Glon[i],Glat[i])
  Rxy=VarT*exp(-Rxy/parest)
  tempre[i]=MT+t(Rxy)%*%V0
}
plotpoints(Glon,Glat,tempre)
lines(Fronti$x,Fronti$y)

############ Tracé dune estimée non paramétrique ############   
########### de la fonction covariance=f(distance) ############
tr=covariog(Dist,T0)
plot(tr$d,tr$r,xlab="Distance",ylab="Covariance")
lines(tr$d,tr$r)
r0 = VarT*exp(-tr$d/parest); #Covariances prédites par le modèle
points(tr$d,r0,pch="x")
lines(tr$d,r0)
title(main="Covariance en fonction de la distance")
legend("topright",pch=c("x","o"),
       legend=c("Prédictions du modèle estimé","Estimées empiriques"))

#########   VALIDATION CROISEE   ##############
######### Calcul de l'erreur par leave one out ##############
Pr=D$data;
cat(" Iteration:")
for (j in 1:NT) {
  cat("",j)
  Dj=D[-j,];
  Dist = 
  VarT=
  MT=
  T0=
  cost = function(par) {

    
    return(ll) 
  }
  par0 #Initialiser intelligemment
  paropt=
  parest=
  R= 
  V0=
  # Calcul de la prediction m+Rxy.V0 en j:
  Rxy=
  Pr[j]=
}
cat("Erreur moyenne: ",sqrt(sum((Pr-D$data)^2)/150))

######### Calcul de l'erreur par leave one out: 
######### algorithme où l'on estime les 3 paramètres au MV

Pr2=D$data;
cat(" Iteration:")
for (j in 1:NT) {
  cat("",j)
  Dj=D[-j,];
  Dist =
  VarT=
  MT=
  T0=
  cost = function(par) {
    R = par[2]*exp(-Dist/par[3]);
    T0=Dj$data-par[1]
    ll = 
    return(ll) 
  }
  par0 = c(xxx,xxx,xxx); Initialiser intelligemment
  paropt=
  parest=
  R=
  T0=
  V0=
  # Calcul de la prediction m+Rxy.V0 en j:
  Rxy=
  Pr2[j]=
}
cat("Erreur moyenne: ",sqrt(sum((Pr2-D$data)^2)/150))


######### BOOTSTRAP par simulation ######### 

library(mvtnorm)

# Caclul du modele 
VarT=var(D$data);
MT=mean(D$data);
T0=D$data-MT;
Dist = lldist(D$lon,D$lat,D$lon,D$lat);
cost = function(par) {
  R = VarT*exp(-Dist/par);
  ll = T0%*%solve(R,T0)+determinant(R)$modulus# log-vraisemblance
  return(ll) 
}
par0 = 160;
paropt=optim(par0,cost,gr=NULL,method="L-BFGS-B")
parest=paropt$par;
R = VarT*exp(-Dist/parest);
NB=100 # Nombre de simulation de bootstrap
parestb=rep(0,NB); # Le vecteur des estimees
Db=D; # Initialisation des donnees de bootstrap
for (j in 1:NB) {
  cat(" ",j)
  Db$data=as.vector(rmvnorm(1,mean=,sigma=)); # Completer
  VarTb=
  MTb=
  T0b=
  cost = function(par) {
 
    return(ll) 
  }
  par0 = 160;
  paropt=optim(par0,cost,gr=NULL,method="L-BFGS-B")
  parestb[j]=paropt$par;
}
hist((parestb-parest)/parest)