[R] Please need help to finalize my code

[PIKAL Petr]

Hi.

Maybe you will get better answers, but from your code it seems to me that you are treating R as C or other similar language which is not optimal. Considering your first 9 lines, it could be changed either to

CenteringV <- function(X, Ms, n) X-Ms

or you probably could use functions sweep or scale.
Your other code is beyond my experience, sorry.

Cheers
Petr

> -----Original Message-----
> From: R-help <r-help-bounces using r-project.org> On Behalf Of Ablaye Ngalaba
> Sent: Saturday, October 10, 2020 9:24 PM
> To: r-help using r-project.org
> Subject: [R] Please need help to finalize my code
> 
> 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 using 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.