#----------------------------------------------------------
#Bibliothèque de fonctions pour le lissage
#----------------------------------------------------------
#Chargement des librairies
library(lattice)
library(latticeExtra)
library(ggplot2)
library(grid)
library(locfit)
theme_set(theme_bw()) 
trellis.device(color = FALSE)
#------------------------------------------------------------------------
resid=function(y,y.fit,infl){
  ##########################
  #Estime les résidus, les résidus normalisés et les résidus de student pour un modèle linéaire
  #'infl' correspond au vecteur de la diagonale de matrice chapeau
  ##########################
  n=length(y)
  res=y-y.fit
  rss=sum(res^2)
  p=sum(infl)
  var.res=rss/(n-p)
  var.res.stud=((n-p)*var.res-res^2/(1-infl))/(n-p-1)
  res.stand=res/sqrt(var.res*(1-infl))
  res.stud=res/sqrt(var.res.stud*(1-infl))
  return(list(resid=res,resid.std=res.stand,resid.stud=res.stud))
}
#------------------------------------------------------------------------
criteria=function(y,y.fit,infl){
  ##########################
  #Estime les différents critères peremettant de déterminer le lissage non paramétrique optimal
  #'infl' correspond au vecteur de la diagonale de matrice chapeau
  ##########################
  n=length(y)
  res=y-y.fit
  rss=sum(res^2)
  p=sum(infl)
  cv=1/n*sum((res/(1-infl))^2)
  gcv=rss/(n*(1-p/n)^2)
  aic=n*log(rss)+2*p
  #bic=n*log(rss)+p*log(n)
  aicc=aic+2*p*(p+1)/(n-p-1)
  return(list(cv=cv,gcv=gcv,aic=aic,aicc=aicc))
}
#------------------------------------------------------------------------
preselect_W_H1=function(y,poids,c_ordre=1:10,c_alpha=1:20){
  ##########################
  #Fonction qui permet par analyse graphique de déterminer les meilleurs paramètres
  #de Whittaker-Henderson en dimension 1 pour l'ensemble des critères de la fonction 'criteria'
  ##########################  
  crit_names=c("cv","gcv","aic","aicc")
  res=expand.grid(crit_names,c_alpha,c_ordre)
  names(res)=c("crit_names","alpha","ordre")
  res$crit=rep(0,nrow(res))
  res$df==rep(0,nrow(res))
  res$crit_names=as.factor(res$crit_names)
  #Fonction critere
  f=function(x){
    y.fit=W_H1(y,poids,x[1],x[2])
    crit=criteria(y,y.fit$fit,y.fit$infl)
    return(c(crit,list(df=y.fit$df)))}
  #alimentation dataset pour graphique
  for(i in c_ordre){
    for(j in c_alpha){
      temp=f(c(i,j))
        for (k in 1:length(crit_names)){
          res[res$crit_names==crit_names[k] & res$alpha==j & res$ordre==i,"crit"]=temp[crit_names[k]]
          res[res$crit_names==crit_names[k] & res$alpha==j & res$ordre==i,"df"]=temp["df"]
        }
    }
  }
  #graphique d'analyse
  graph=ggplot(res,aes(x=alpha,crit,label=ordre))+
    geom_text(size=3)
  graph=graph+facet_wrap(~crit_names,2,2,scales="free_y")
  print(graph)
}
#------------------------------------------------------------------------
compar_W_H1=function(x,y,poids,c_ordre,c_alpha){
  ##########################
  #Fonction qui permet par analyse graphique de déterminer les paramètres optimaux
  #de Whittaker-Henderson en dimension 1 pourl'ensemble des critères de la fonction 'criteria'
  ##########################
  res=expand.grid(y=y,ordre=c_ordre,alpha=c_alpha)
  res$y.fit=0
  res$resid=0
  res$resid.fit=0
  res$x=rep(x,length(c_ordre)*length(c_alpha))
  for(i in c_ordre){
    for(j in c_alpha){
      y.fit=W_H1(y,poids,i,j)
      res.brut=resid(y,y.fit$fit,y.fit$infl)$resid
      res.fit=locfit(res.brut~x, family = "gaussian", alpha = 0.2, link = "identity", kern="gauss")    
      res$y.fit[res$ordre==i & res$alpha==j]=y.fit$fit
      res$resid[res$ordre==i & res$alpha==j]=res.brut
      res$resid.fit[res$ordre==i & res$alpha==j]=fitted(res.fit)
      rm(y.fit,res.brut,res.fit)
    }
  }
  #graphiques d'analyses
  #valeurs lissées
  graph=ggplot(res,aes(x=x,y=y))+
    geom_point(shape=1)+
    geom_line(aes(x=x,y=y.fit))
  graph=graph+
    scale_x_continuous("Age")+
    scale_y_continuous("Qx")+
    facet_grid(alpha ~ ordre,scales="free_y")
  print(graph)
  #résidus
  graph=ggplot(res,aes(x=x,y=resid))+
    geom_point(shape=1)+
    geom_line(aes(x=x,y=resid.fit),size=1)+
    geom_hline(yintercept=0,linetype=2)
  graph=graph+
    scale_x_continuous("Age")+
    scale_y_continuous("Residus")+
    facet_grid(alpha ~ ordre,scales="free_y")
  print(graph)
}
#------------------------------------------------------------------------
plot.res.liss=function(data,poids,names_trans){
  ##########################
  #Fonction qui de tracer les résultats du lissage
  #Les taux bruts et lissés sont comparés ainsi que les résidus
  #Chaque jeu de données est codé par la colonne 'trans'du data.frame
  #Le data.frame doit avoir les colonnes "Age","qx.brut","trans","qx.fit","infl","resid","resid.std","resid.stud"
  ##########################
  n=length(names_trans)
  for(i in 1:n){
    temp=data[data$trans==names_trans[i],]
    res.fit=locfit(temp$resid~temp$Age, family = "gaussian", alpha = 0.2, link = "identity", kern="gauss")    
    data$resid.fit[data$trans==names_trans[i]]=fitted(res.fit)
    rm(temp,res.fit)
    }
  #Graphiques
  plot_1=ggplot(data,aes(x=Age,y=qx.brut))+
    geom_point(shape=1)+
    geom_line(aes(x=Age,y=qx.fit))+
    scale_x_continuous("Age")+
    scale_y_continuous("Taux d'incidence")+
    facet_wrap(~trans,ncol=1,scales="free_y")
  plot_2=ggplot(data,aes(x=Age,y=resid))+
    geom_point(shape=1)+
    geom_line(aes(x=Age,y=resid.fit),size=1)+
    scale_x_continuous("Age")+
    scale_y_continuous("Residus")+
    facet_wrap(~trans,ncol=1,scales="free_y")
  plot_3=ggplot(data,aes(x=Age,y=resid.stud))+
    geom_point(shape=1)+
    geom_hline(yintercept=0, linetype=2)+
    geom_hline(yintercept=2, linetype=3)+
    geom_hline(yintercept=-2, linetype=3)+
    scale_x_continuous("Age")+
    scale_y_continuous("Residus Student",limits=c(-4,4))+
    facet_wrap(~trans,ncol=1,scales="free_x")
  return(grid.arrange(plot_1,plot_2,plot_3,nrow=1,ncol=3))
}
#------------------------------------------------------------------------
test_process.fit=function(data,names_trans,t.stat){
  ##########################
  #Fonction qui permet d'exécuter une batterie de tests pour analyser la qualité des ajustements
  #Fonctionne sur un dataframe de résultats pour lequel chaque jeu de données est codé par 
  #la colonne 'trans'du data.frame
  #Le data.frame doit avoir les colonnes "Age","qx.brut","trans","qx.fit","infl","resid","resid.std","resid.stud"
  ##########################
  #1.Degré de liberté
  print("### Degre de liberte ###")
  res.deg=matrix(,1,length(names_trans))
  colnames(res.deg)=names_trans
  rownames(res.deg)=c("Degre de liberte")
  for (i in 1:length(names_trans)){
    res.deg[1,i]=sum(data$infl[data$trans==names_trans[i]])
  }
  print(res.deg)
  #2.Residus stud
  print ("### Residus de student ###")
  res.resid.stud=matrix(,2,length(names_trans))
  colnames(res.resid.stud)=names_trans
  rownames(res.resid.stud)=c(">2",">3")
  for (i in 1:length(names_trans)){
    res.resid.stud[1,i]=length(which(abs(data$resid.stud[data$trans==names_trans[i]])>2))
    res.resid.stud[2,i]=length(which(abs(data$resid.stud[data$trans==names_trans[i]])>3))
  }
  print(res.resid.stud)
  #3.Kolmogorov-Smirnov
  print ("### Komogorov-Smirnov ###")
  res.ks=matrix(,2,length(names_trans))
  colnames(res.ks)=names_trans
  rownames(res.ks)=c("Statistique","p.value")
  for (i in 1:length(names_trans)){
    ks=ks.test(data$qx.brut[data$trans==names_trans[i]],data$qx.fit[data$trans==names_trans[i]])
    res.ks[1,i]=ks$statistic
    res.ks[2,i]=ks$p.value
    rm(ks)
  }
  print(res.ks)
  #4.Chi2
  print ("### Chi 2 ###")
  res.chi2=matrix(,2,length(names_trans))
  colnames(res.chi2)=names_trans
  rownames(res.chi2)=c("Statistique","p.value")
  for (i in 1:length(names_trans)){
    expo=t.stat$expo[t.stat$Trans==names_trans[i]]
    y=data$qx.brut[data$trans==names_trans[i]]
    y.fit=data$qx.fit[data$trans==names_trans[i]]
    chi=sum(((y.fit-y)^2)*expo/(y*(1-y)))
    res.chi2[1,i]=chi
    res.chi2[2,i]=pchisq(chi,df=length(y.fit)-res.deg[1,i]-1)
    rm(expo,y,y.fit,chi)
  }
  print(res.chi2)
  #5. R^2
  print ("### R^2 ###")
  res.R2=matrix(,1,length(names_trans))
  colnames(res.R2)=names_trans
  rownames(res.R2)=c("valeur")
  for (i in 1:length(names_trans)){
    y=data$qx.brut[data$trans==names_trans[i]]
    e=data$resid[data$trans==names_trans[i]]
    num=sum(e^2)
    denum=sum((y-mean(y))^2)
    res.R2[1,i]=1-num/denum
    rm(y,e,num,denum)
  }
  print(res.R2)
  return(list(deg=res.deg,resid.stud=res.resid.stud,ks=res.ks,chi2=res.chi2,R2=res.R2))
}
#-----------------------------------------------------------------------------------------
#Lissage par Whittaker-Henderson en dimension 1
W_H1<-function(y,Poids,Ordre,alpha){ 
  #Déclaration de la taille des variables
  q=length(y)
  #Définition des matrices
  U=matrix(0,q)
  v=matrix(0,Ordre+1)
  Kv=matrix(0,(q-Ordre),q)
  M=matrix(0,q,q)
  fit=matrix(0,q)
  #Construction de la matrice des poids
  W=diag(Poids)
  #Construction matrice Kv
  for (k in 0:Ordre) {
    v[(k+1),1]=(-1)^(Ordre-k)*factorial(Ordre)/(factorial(k)*factorial(Ordre-k))
  }
  for (z in 1:(q-Ordre)) {
    for (i in 1:(Ordre+1)) { 
      Kv[z,i+(z-1)]=v[i,1]
    } 
  } 
  #Calcul des taux lissés
  M=W+alpha*t(Kv)%*%Kv
  H=solve(M)%*%W
  fit=H%*%y
  infl=diag(H)
  df=sum(infl)
  return(list(fit=fit,infl=infl,df=df))
}
#-----------------------------------------------------------------------------------------
#Lissage par Whittaker-Henderson en dimension 2
W_H2<-function(Tbrut,Poids,Ordre_vertical,Ordre_horizontal,alpha,beta){
  
  #Déclaration de la taille des variables
  q<<-ncol(Tbrut)
  p<<-nrow(Tbrut)
  
  #Définition des matrices
  U=matrix(0,q*p)
  v=matrix(0,Ordre_vertical+1)
  h=matrix(0,Ordre_horizontal+1)
  x=matrix(1:p,p,1)
  y=matrix(1:q,q,1)
  Kv=matrix(0,(p-Ordre_vertical)*q,q*p)
  Kh=matrix(0,p*(q-Ordre_horizontal),q*p)
  M=matrix(0,q*p,q*p)
  W=matrix(0,q*p,q*p)
  qlisse=matrix(0,p,q)
  Tlisse=matrix(0,p,q)
  Ecart=matrix(0,p,q)
  
  #Transformation de la matrice de taux bruts en vecteur et construction de la matrice des poids
  for (j in 1:q) {
    for (i in 1:p) {
      U[(i-1)*q+j,1]=Tbrut[i,j]
      W[q*(i-1)+j,q*(i-1)+j]=Poids[i,j]
    } 
  }
  
  #Construction matrice kv
  for (k in 0:Ordre_vertical) {
    v[(k+1),1]=(-1)^(Ordre_vertical-k)*factorial(Ordre_vertical)/(factorial(k)*factorial(Ordre_vertical-k))
  }
  
  for (j in 1:q) {
    for (z in 1:(p-Ordre_vertical)) {
      for (i in 1:(Ordre_vertical+1)) { 
        Kv[z+(j-1)*(p-Ordre_vertical),j+(q)*(i-1)+(z-1)*(q)]=v[i,1]
      } 
    } 
  }
  
  #Construction matrice Kh
  for (k in 0:Ordre_horizontal) {
    h[(k+1),1]=(-1)^(Ordre_horizontal-k)*factorial(Ordre_horizontal)/(factorial(Ordre_horizontal-k)*factorial(k))
  }
  for (i in 1:p) {
    for (j in 1:(q-Ordre_horizontal)) {
      for (z in 1:(Ordre_horizontal+1)) { 
        Kh[j+(i-1)*(q-Ordre_horizontal),z+(j-1)+(i-1)*(q)]=h[z,1]
      }
    }
  }
  
  #Calcul des taux lissés
  M=W+alpha*t(Kv)%*%Kv+beta*t(Kh)%*%Kh
  qlisse=solve(M)%*%W%*%U
  for (j in 1:q) {
    for (i in 1:p) {
      Tlisse[i,j]=qlisse[(i-1)*q+j,1]
      Ecart[i,j]=Tlisse[i,j]/Tbrut[i,j]
    } 
  }
  
  #Stockage des résultats
  Tlisse<-Tlisse
  Ecart<-Ecart
  
  #Représentation graphique
  par(mfrow=c(1,3))
  persp(x,y,Tbrut,theta=45,col=brewer.pal(9,"RdYlGn"),xlab="Age",ylab="Ancienneté",main="Taux bruts",adj=0.5,font=2)
  persp(x,y,Tlisse,theta=45,col=brewer.pal(9,"RdYlGn"),xlab="Age",ylab="Ancienneté",main="Taux lissés",adj=0.5,font=2)
  persp(x,y,Ecart,theta=90,phi=22.5,expand=1,col=brewer.pal(9,"RdYlGn"),xlab="Age",ylab="Ancienneté",main="Taux lissés/Taux bruts",adj=0.5,font=2)
  
}