[R] Aide pour finaliser ce code

[Ablaye Ngalaba]

Hello.
Here is my R code. I used the functional data . Now I need to use the
functional data by applying the kernels instead of the xi, yi functions.


Bonjour.
Voici mon code en R . J'ai utiliser les données fonctionnelles . Maintenant
j'ai besoin d'utiliser les données fonctionnelles en appliquant les noyaux
à la place des fontions xi, yi





library(MASS)

CentrageV<-function(X,Ms,n){
# cette fonction centre les données de X
X1=X*0
for (i in 1:n){
X1[i,]=X[i,]-Ms
}
return(X1)
}


# Fonction N°2
SqrtMatrice0<-function(M) {
# Cette fonction nous permet de calculer la racine carrée d'une matrice
# en utilisant la décomposition M=PDQ où Q est l'inverse de P
# ici les valeurs propres négatives sont remplacées par zero
 a=eigen(M,TRUE)
 b=a$values
 b[b<0]=0
 c=a$vectors
 d=diag(sqrt(b))
 b=solve(c)
 a=c%*%d%*%b
return(a)
}


# déclaration des parametres
m1=0.01 # valeur de alpha (risque de 1%)
m2=0.05 # valeur de alpha (risque de 5%)
m3=0.1 # valeur de alpha (risque de 10%)
nbrefoissim=100 # nbrefois que le programme tourne
p=2 #dimension de la variable X
q=3 #dimension de la variable Y
R=c(2,3,2);# Nbre de partition de chaque composante de la variable Y
if(length(R) != q) stop("La taille de R doit être égale à q")
n=25 # Taille echantillon
N=c(25,50,100,200,300,400,500,1000) #differentes tailles de l'échantillion
N=c(25,50) #differentes tailles de l'échantillion

K=0
MV=matrix(0,nr=2,nc=4)
dimnames(MV)[[2]]=c("n","r1%","r5%","r10%")


#Debut du programme

 for (n in N){


l1=0 # initialisation de la valeur permettant calculer le niveau de test à
1%
l2=0 # initialisation de la valeur permettant calculer le niveau de test à
5%
l3=0 # initialisation de la valeur permettant calculer le niveau de test à
10%

# Création d'une liste n11 qui contient les tailles des differents groupes
n11=list()
for (i in 1:q){
n11[[i]]=rep(as.integer(n/R[i]),R[i])
n11[[i]][R[i]]=n-((R[i]-1)*n11[[i]][1])
}


# Création des listes  P11 et P12 qui contient les probabilités et
# les inverses des probabilites empiriques des differents groupes
respectivement

P11=list()
P12=list()
for (i in 1:q){
P11[[i]]=n11[[i]]/n
P12[[i]]=n/n11[[i]]

}

# création d'une liste contenant les matrices W
W=list()
for (i in 1:q){
w=matrix(0,n,R[i])
w[1:n11[[i]][1],1]=1
for (j in 2:R[i]){
s1=sum(n11[[i]][1:(j-1)])
w[(1+s1):(s1+n11[[i]][j]),j]=1
}
W[[i]]=w
}

for (i1 in 1:nbrefoissim){

# géneration des données
VA1=mvrnorm(n,rep(0,(p+q)),diag((p+q)))
X=VA1[,1:p]
Y=VA1[,(p+1):(p+q)]

# Calcul de Xbar
Xbar=colMeans(X)

# Calcul des Xjh bar
Xjhbar=list()
for (i in 1:q){
  w=matrix(0,R[i],p)
  for (j in 1:R[i]){
     w[j,]=colSums(W[[i]][,j]*X)/n11[[i]][j]
  }
  Xjhbar[[i]]=w
}

#calcul des TO jh

TO.jh=list()
for (i in 1:q){
  w=Xjhbar[[i]]
  to=w*0
  for (j in 1:R[i]){
     to[j,]=w[j,]-Xbar
  }
  TO.jh[[i]]=to
}

#calcul des Lamda J

Lamda=matrix(0,p,p)
for (i in 1:q){
  to=TO.jh[[i]]
  w=matrix(0,p,p)
  for (j in 1:R[i]){
     w=w+(P11[[i]][j]*(to[j,]%*%t(to[j,])))
  }
  Lamda=Lamda+w
}

tr1=n*sum(diag(Lamda))

# Calcul de Gamma

GGamma=matrix(0,p*sum(R),p*sum(R))
PGamma=kronecker(diag(P11[[1]]),diag(p))
Ifin=p*R[1]
GGamma[1:Ifin,1:Ifin]=PGamma
for (i in 2:q){
  PGamma=kronecker(diag(P11[[i]]),diag(p))
  Idebut=((p*sum(R[1:(i-1)]))+1)
  Ifin=(p*sum(R[1:i]))
  GGamma[Idebut:Ifin,Idebut:Ifin]=PGamma
}

#Calcul de Sigma

  # Calcul de Vn
  X1=CentrageV(X,Xbar,n)
  Vn=t(X1)%*%X1/n

## Construction de Sigma
GSigma=matrix(0,p*sum(R),p*sum(R))
for (i in 1:q ){
   for (j in 1:R[i] ){
      for (k in 1:q ){
         for (l in 1:R[k]){
            Xij=CentrageV(X,Xjhbar[[i]][j,],n)
            Xkl=CentrageV(X,Xjhbar[[k]][l,],n)
            Vijkl=t(W[[i]][j]*Xij)%*%(W[[k]][l]*Xkl)/n
            Vij=t(W[[i]][j]*Xij)%*%Xij/n11[[i]][j]
            Vkl=t(W[[k]][l]*Xkl)%*%Xkl/n11[[k]][l]
            if (i==1) Idebut=((j-1)*p)+1 else
Idebut=((sum(R[1:(i-1)])+j-1)*p)+1
   if (i==1) Idebut1=(j*p) else Idebut1=((sum(R[1:(i-1)])+j)*p)
            if (k==1) Ifin=((l-1)*p)+1 else Ifin=((sum(R[1:(k-1)])+l-1)*p)+1
   if (k==1) Ifin1=(l*p) else Ifin1=((sum(R[1:(k-1)])+l)*p)

 GSigma[Idebut:Idebut1,Ifin:Ifin1]=Vijkl+(P11[[i]][j]*P11[[k]][l]*(Vn-Vij-Vkl))
}
     }
   }
 }


# Déterminations des valeurs propres de sigma epsilon
pa=SqrtMatrice0(GSigma)
mq= pa %*% GGamma %*% pa
u=Re(eigen(mq)$values)


# détermination de dégré de liberté et valeur c noté va
dl=(sum(u)^2)/(sum(u^2))
va=(sum(u^2))/(sum(u))
pc=1-pchisq(tr1/va, df= dl)

# Test de la valeur obtenue
if (pc>m1) d1=0 else d1=1
if (pc>m2) d2=0 else d2=1
if (pc>m3) d3=0 else d3=1
l1=l1+d1
l2=l2+d2
l3=l3+d3
}
K=K+1
MV[K,1]=n
MV[K,2]=l1/nbrefoissim
MV[K,3]=l2/nbrefoissim
MV[K,4]=l3/nbrefoissim
}

	[[alternative HTML version deleted]]