Aide pour finaliser ce code

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

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]]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Aide pour finaliser ce code

Jeremie Juste
Hello Ablaye,

I didn't find any function xi ,yi.

This is an english mailing list. You will maximize your chance of
getting help by sticking to english.

I would also recommend that you comment your code in English, it will be
easier to share. But on this one it is your call.

PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Best regards,
Jeremie Juste

On Friday,  9 Oct 2020 at 13:31, Ablaye Ngalaba wrote:

> 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]]
>
> ______________________________________________
> [hidden email] 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.
>

______________________________________________
[hidden email] 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.