[R] Please need help to finalize my code
PIKAL Petr
petr@p|k@| @end|ng |rom prechez@@cz
Wed Oct 14 10:35:24 CEST 2020
Hi
I do not want to use kernels, OP wants, God knows for what.
Cheers
Petr
From: Richard O'Keefe <raoknz using gmail.com>
Sent: Wednesday, October 14, 2020 1:43 AM
To: PIKAL Petr <petr.pikal using precheza.cz>
Cc: Ablaye Ngalaba <ablayengalaba using gmail.com>; r-help using r-project.org
Subject: Re: [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 http://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 <mailto:petr.pikal using 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 <mailto:ablayengalaba using gmail.com>
Sent: Monday, October 12, 2020 7:58 PM
To: PIKAL Petr <mailto:petr.pikal using 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 <mailto:petr.pikal using precheza.cz <mailto:mailto:petr.pikal using 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 <mailto:r-help-bounces using r-project.org> On Behalf Of Ablaye Ngalaba
> Sent: Saturday, October 10, 2020 9:24 PM
> To: mailto:r-help using r-project.org <mailto:mailto: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]]
>
> ______________________________________________
> mailto:R-help using r-project.org <mailto:mailto: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.
______________________________________________
mailto: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.
More information about the R-help
mailing list