Good evening dear administrators,
It is with pleasure that I am writing to you to ask for help to finalize my
R programming algorithm.
Indeed, I attach this note to my code which deals with a case of
independence test statistic . My request is to introduce the kernels using
the functional data for this same code that I am sending you. So I list the
lines for which we need to introduce the functional data using the kernels
below.
First of all for this code we need to use the functional data. I have
numbered the lines myself from the original code to give some indication
about the lines of code to introduce the kernel.
1)Centering of redescription function phi(xi) (just use the kernel)
2)this function centers the functional data phi(xi) (just use the kernel)
3)phi(x1) (or kernel k( ., .) )
5)Even here the kernel with the functional data is used.
7)return the kernel
28) kernel dimension
29) kernel dimension
30)kernel dimension
73) redescription function phi(xi) or the kernel k(. , .)
74)redescription function phi(yi) or the kernel k(. , .)
76)redescription function phi(xbar) or the kernel k(. , .)
77)redescription function phi(xjhbar) or the kernel k(. , .)
82)redescription function phi(xi) or the kernel k(. , .) introduced instead
of x
84)redescription function phi(xjhbar) or the kernel k(. , .)
89)w= redescription function phi(xbar) or the kernel k(. , .)
120) we center the redescription function phi(xi) or the kernel k(. , .)
121)Vn is a kernel function, so we use ph(xi) or the kernel k(. , , .)
126)centering of the redescription function phi(xji) or the kernel k(. , .)
127)Vij is computed using redescription function phi(xi) or the kernel k(.
, .) instead of X.
130)centering redescription function phi(Xkl) or the kernel k(. , .)
131)Vijkl always using the kernel function as precedent
132)Vkl always using the kernel function as precedent
I look forward to your help and sincere request.
Yours sincerely
library(MASS)
1 CenteringV<-function(X,Ms,n){
2 # this function centers the data of X
3 X1=X*0
4 for (i in 1:n){
5 X1[i,]=X[i,]-Ms
6 }
7 return(X1)
8 }
9 # Function N°2
10 SqrtMatrix0<-function(M) {
11 # This function allows us to compute the square root of a matrix
12 # using the decomposition M=PDQ where Q is the inverse of P
13 # here negative eigenvalues are replaced by zero
14 a=eigen(M,TRUE)
15 b=a$values
16 b[b<0]=0
17 c=a$vectors
18 d=diag(sqrt(b))
19 b=solve(c)
20 a=c%*%d%*%b
21 return(a)
22 }
23 # declaration of parameters
24 m1=0.01 # alpha value (1% risk)
25 m2=0.05 # alpha value (5% risk)
26 m3=0.1 # alpha value (10% risk)
27 nbrefoissim=100 # # times the program is running
28 p=3 #dimension of variable X
29 q=3 #dimension of variable Y
30 R=c(5,4,3);# Number of partition of each component of variable Y
31 if(length(R) != q) stop("The size of R must be equal to q")
32 n=25 # Sample size
33 N=c(25,50,100,200,300,400,500,1000) #different sample sizes
34 #N=c(25,50) #different sample sizes
35 pc1=rep(0.100)
36 K=0
37 MV=matrix(0,nr=8,nc=4)
38 dimnames(MV)[[2]]=c("n","r1%","r5%","r10%")
39 #Beginning of the program
40 for (n in N){
41 l1=0 # initialization of the value to calculate the test level at 1%.
42 l2=0 # initialization of the value to calculate the test level at 5%.
43 l3=0 # initialization of the value to calculate the test level at 10%.
44 # Creation of an n11 list containing the sizes of the different groups
45 n11=list()
46 for (i in 1:q){
47 n11[[i]]=rep(as.integer(n/R[i]),R[i])
48 n11[[i]][R[i]]=n-((R[i]-1)*n11[[i]][1])
49 }
50 # Creation of lists P11 and P12 which contain the probabilities and
51 # the inverses of the empirical probabilities of the different groups
respectively
52 P11=list()
53 P12=list()
54 for (i in 1:q){
55 P11[[i]]=n11[[i]]/n
56 P12[[i]]=n/n11[[i]
57 }
58 # creation of a list containing the W matrices
59 W=list()
60 for (i in 1:q){
61 w=matrix(0,n,R[i])
62 w[1:n11[[i]][1],1]=1
63 if (R[i]>1){
64 for (j in 2:R[i]){
65 s1=sum(n11[[i]][1:(j-1)])
66 w[(1+s1):(s1+n11[[i]][j]),j]=1
67 }}
68 W[[i]]=w
69 }
70 for (i1 in 1:nbrefoissim){
71 # data generation
72 VA1=mvrnorm(n,rep(0,(p+q)),diag((p+q)))
73 X=VA1[,1:p]
74 Y=VA1[,(p+1):(p+q)]
75 # Xbar calculation
76 Xbar=colMeans(X)
77 # Calculation of Xjh bar
78 Xjhbar=list()
79 for (i in 1:q){
80 w=matrix(0,R[i],p)
81 for (j in 1:R[i]){
82 w[j,]=colSums(W[[i]][,j]*X)/n11[[i]][j]
83 }
84 Xjhbar[[i]]=w
85 }
86 #calculation of TO jh
87 TO.jh=list()
88 for (i in 1:q){
89 w=Xjhbar[[i]]
90 to=w*0
91 for (j in 1:R[i]){
92 to[j,]=w[j,]-Xbar
93 }
94 TO.jh[[i]]=to
95 }
96 #calculation of Lamda J
97 Lamda=matrix(0,p,p)
98 for (i in 1:q){
99 to=TO.jh[[i]]
100 w=matrix(0,p,p)
101 for (j in 1:R[i]){
102 w=w+(P11[[i]][j]*(to[j,]%*%t(to[j,])))
103 }
104 Lamda=Lamda+w
105 }
106 tr1=n*sum(diag(Lamda))
107 # Gamma Calculation
108 GGamma=matrix(0,p*sum(R),p*sum(R))
109 PGamma=kronecker(diag(P12[[1]]),diag(p))
110 Ifin=p*R [1]
111 GGamma[1:Ifin,1:Ifin]=PGamma
112 for (i in 2:q){
113 PGamma=kronecker(diag(P12[[i]]),diag(p))
114 Idebut=((p*sum(R[1:(i-1)]))+1)
115 Ifin=(p*sum(R[1:i]))
116 GGamma [Idebut:Ifin,Idebut:Ifin]=PGamma
117 }
118 #Sigma calculation
119 # Calculation of Vn
120 X1=CenteringV(X,Xbar,n)
121 Vn=t(X1)%*%X1/n
122 ## Construction of Sigma
123 GSigma=matrix(0,p*sum(R),p*sum(R))
124 for (i in 1:q ){
125 for (j in 1:R[i] ){
126 Xij=CenteringV(X,Xjhbar[[i]][j,],n)
127 Vij=t(W[[i]][,j]*Xij)%*%Xij/n11[[i]][j]
128 for (k in 1:q ){
129 for (l in 1:R[k]){
130 Xkl=CenteringV(X,Xjhbar[[k]][l,],n)
131 Vijkl=t(W[[i]][,j]*Xij)%*%(W[[k]][,l]*Xkl)/n
132 Vkl=t(W[[k]][,l]*Xkl)%*%Xkl/n11[[k]][l]
133 if (i==1) Idebut=((j-1)*p)+1 else Idebut=((sum(R[1:(i-1)])+j-1)*p)+1
134 if (i==1) Idebut1=(j*p) else Idebut1=((sum(R[1:(i-1)])+j)*p)
135 if (k==1) Ifin=((l-1)*p)+1 else Ifin=((sum(R[1:(k-1)])+l-1)*p)+1
136 if (k==1) Ifin1=(l*p) else Ifin1=((sum(R[1:(k-1)])+l)*p)
137
GSigma[Idebut:Idebut1,Ifin:Ifin1]=Vijkl+(P11[[i]][j]*P11[[k]][l]*(Vn-Vij-Vkl))
138 }
139 }
140 }
141 }
141 # Determinations of the eigenvalues of sigma epsilon
142 pa=SqrtMatrix0(GSigma)
143 mq= pa %*% GGamma %*% pa
144 u=Re(eigen(mq)$values)
145 # determination of degree of freedom and value c noted va
146 dl=(sum(u)^2)/(sum(u^2))
147 va=(sum(u^2))/(sum(u))
148 pc=1-pchisq(tr1/va, df= dl)
149 pc1[i1]=pc
150 # Test of the value obtained
151 if (pc>m1) d1=0 else d1=1
152 if (pc>m2) d2=0 else d2=1
153 if (pc>m3) d3=0 else d3=1
154 l1=l1+d1
155 l2=l2+d2
156 l3=l3+d3
157 }
158 K=K+1
159 MV [K,1]=n
160 MV [K,2]=l1/number of timessim
161 MV [K,3]=l2/number of timessim
162 MV[K,4]=l3/number of timessim
163 }
l
[[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.
