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]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
