thr3ads.net - R help - [R] Please need help to finalize my code [Oct 2020]

If this information is useful, please help other people find it:
Share via:

Ablaye Ngalaba

2020-Oct-10 19:24 UTC

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

PIKAL Petr

2020-Oct-12 06:28 UTC

head link

[R] Please need help to finalize my code

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 at r-project.org> On Behalf Of Ablaye
Ngalaba
> Sent: Saturday, October 10, 2020 9:24 PM
> To: r-help at 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 at 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.

PIKAL Petr

2020-Oct-13 05:21 UTC

head link

[R] Please need help to finalize my code

Hm. Google tells me that kernel function is in stats package which comes with
base installation and is invoked when you start R.

 

search()

[1] ".GlobalEnv"        "package:stats"    
"package:graphics"

[4] "package:grDevices" "package:utils"    
"package:datasets"

[7] "package:methods"   "Autoloads"        
"package:base"

 

And BTW, keep your posts on R-help, others could give you more relevant answers.

 

Cheers

Petr

 

From: Ablaye Ngalaba <ablayengalaba at gmail.com> 
Sent: Monday, October 12, 2020 7:58 PM
To: PIKAL Petr <petr.pikal at precheza.cz>
Subject: Re: [R] Please need help to finalize my code

 

Hello.
Thank you for your response.
I would like to ask you the name of the package to install when you want to use
the kernels

 

Le lun. 12 oct. 2020 ? 08:28, PIKAL Petr <petr.pikal at precheza.cz
<mailto:petr.pikal at precheza.cz> > a ?crit :

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 at r-project.org> On Behalf Of Ablaye
Ngalaba
> Sent: Saturday, October 10, 2020 9:24 PM
> To: r-help at r-project.org <mailto:r-help at 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 at r-project.org <mailto:R-help at 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.

Richard O'Keefe

2020-Oct-13 23:43 UTC

head link

[R] Please need help to finalize my code

What do you *mean* "when you want to use the kernels".
WHICH kernels?
Use to do WHAT?
In your browser, visit cran.r-project.org
then select "Packages" from the list on the left.
Then pick the alphabetic list.
Now search for 'kernel'.
You will find dozens of matches.

On Wed, 14 Oct 2020 at 05:15, PIKAL Petr <petr.pikal at precheza.cz>
wrote:
> Hm. Google tells me that kernel function is in stats package which comes
> with base installation and is invoked when you start R.
>
>
>
> search()
>
> [1] ".GlobalEnv"        "package:stats"    
"package:graphics"
>
> [4] "package:grDevices" "package:utils"    
"package:datasets"
>
> [7] "package:methods"   "Autoloads"        
"package:base"
>
>
>
> And BTW, keep your posts on R-help, others could give you more relevant
> answers.
>
>
>
> Cheers
>
> Petr
>
>
>
> From: Ablaye Ngalaba <ablayengalaba at gmail.com>
> Sent: Monday, October 12, 2020 7:58 PM
> To: PIKAL Petr <petr.pikal at precheza.cz>
> Subject: Re: [R] Please need help to finalize my code
>
>
>
> Hello.
> Thank you for your response.
> I would like to ask you the name of the package to install when you want
> to use the kernels
>
>
>
> Le lun. 12 oct. 2020 ? 08:28, PIKAL Petr <petr.pikal at precheza.cz
<mailto:
> petr.pikal at precheza.cz> > a ?crit :
>
> 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 at r-project.org> On Behalf Of
Ablaye Ngalaba
> > Sent: Saturday, October 10, 2020 9:24 PM
> > To: r-help at r-project.org <mailto:r-help at 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 at r-project.org <mailto:R-help at 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.
>
> ______________________________________________
> R-help at 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.
>
	[[alternative HTML version deleted]]

Reasonably Related Threads

Search for more seemingly similar threads

R help - Oct 2020 - Please need help to finalize my code

[R] Please need help to finalize my code

[R] Please need help to finalize my code

[R] Please need help to finalize my code

[R] Please need help to finalize my code

Reasonably Related Threads