Hello,
I need help with my R programming code.
I just have a problem storing the matrix calculate Gsigma[m] in sigma.
Seeing what I have coded I find NA value which doesn't generate me errors.
library(MASS)
# Creation of the linear kernel function
#Function N°1
Kernellin<-function(X,Y){
return(X%*%t(Y))
}
# Function N°2
SqrtMatrice0<-function(M) {
# This function allows us to calculate the square root of a matrix
# using the decomposition M=PDQ where Q is the inverse of P
# here the negative eigenvalues are replaced by 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)
}
# declaration of parameters
m1=0.01 # value of alpha (risque de 1%)
m2=0.05 # value of alpha (risque de 5%)
m3=0.1 # value of alpha (risque de 10%)
nbrefoissim=100 # number of times the program runs
p=3 #dimension of the variable X
#q=3 #dimension of the variable Y
#R=c(5,4,3);# Number of partitions of each component of the variable Y
if(length(R) != q) stop("The size of R must be equal to q")
n=25 # Taille echantillon
N=c(25,50,100,200,300,400,500,1000) #different sample sizes
#N=c(25,50) #ifferent sample sizes
pc1=rep(0,100)
K=0
MV=matrix(0,nr=8,nc=4)
dimnames(MV)[[2]]=c("n","r1%","r5%","r10%")
#Start of the programme
for (n in N){
l1=0 # initialization of the value to calculate the test level at 1%.
l2=0 # initialization of the value to calculate the test level at 1%. 5%
l3=0 # initialization of the value to calculate the test level at 1%. 10%
# Creation of a list nl which contains the sizes of the different groups
y <- sample(q, 10, TRUE)
l <- sample(q, 10, TRUE)
ind <- function(y,l) ifelse(y == l, 1,0)
for (i in 1:q) {
nl[[i]]<-rep(sum(ind(y,l)),i)
nl[[i]][[i]]=n-((i-1)*nl[[i]][1])
}
# Creation of the lists Pl1 and Pl2 which contain the probabilities and
# the inverses of the empirical probabilities of the different groups
respectively
pl1<-list()
pl2<-list()
for (i in 1:p) {
pl1[[i]]<-nl[[i]]/n
pl2[[i]]<-n/nl[[i]]
}
#Creation of a matrix unit
Tt<-matrix (rep(1, 3*3), 3, 3)
#Création de matrice
aa<-c(1,0,1)
bb<-c(1,1,0)
TT<-matrix(c(aa,bb),3,3)
for (i1 in 1:nbrefoissim){
# data generation
X <- mvrnorm(3,rep(0,p),diag(p))
#Sigma calculation
#kernel gram matrix calculation
K1<-TT%*%Kernellin(X,X)%*%t(TT)
k2<-Tt%*%Kernellin(X,X)%*%t(Tt)
k3<-TT%*%Kernellin(X,X)%*%t(Tt)
k4<-Tt%*%Kernellin(X,X)%*%t(TT)
# B calculation
B <- matrix(0, p, p)
for (l in 1:q) {
B<-B + pl1[[l]][1]*( K1/n^2 - k3/n^2 - k4/n^2 + k2/n^2 )
}
#kernel gram matrix calculation
K5<-t(Tt)%*%Kernellin(X,X)
K6<-Kernellin(X,X)%*%Tt
#V calculation
V <- matrix(0, p, p)
V<-V + (Kernellin(X,X) - K5/n - K6/n - k4/n^2 )/n
#kernel gram matrix calculation
K7<-TT%*%Kernellin(X,X)
K8<-Kernellin(X,X)%*%TT
#W calculation
W<-list()
for (i in 1:p) {
for (j in 1:p) {
W<-matrix(0, p, p)
Wi<-W + (Kernellin(X,X) - K7/n - K8/n - K1/n^2 )/nl[[i]][1]
Wj<-W + (Kernellin(X,X) - K7/n - K8/n - K1/n^2 )/nl[[j]][3]
}
}
# T4 calculation
T4<-n*sum(diag(V%*%B))
# GGamma calculation
GGamma<-matrix(0,p*p,p*p)
GGamma<- kronecker(diag(pl2[[i]]), ginv(V))
#GSigma calculation
kron <- function(i, j) ifelse(i==j, 1, 0)
GSigma <- list()
m = 1
for(i in 1:p){
for(j in 1:p){
GSigma[[m]] <- kron(i,j)*pl1[[l]][1]*Wi +pl1[[i]][1]*pl1[[j]][3]*(V
- Wi - Wj )
m <- m+1
}
}
sigma <- matrix(0, p*p, p*p)
m <- 1
s <- 1
for(i in 1:p){
r <- i*p
a <- 1
for(j in 1:p){
b <- j*p
sigma[s:r, a:b] <- GSigma[[m]]
m <- m+1
a <- b+1
}
s <- r+1
}
# Eigenvalue determinations of sigma epsilon
pa<-SqrtMatrice0(sigma)
mq<- pa %*% GGamma %*% pa
u<-Re(eigen(mq)$values)
# determination of degree of freedom and c-value noted va
dl<-(sum(u)^2)/(sum(u^2))
va<-(sum(u^2))/(sum(u))
pc<-1-pchisq(tr1/va, df= dl)
pc1[i1]<-pc
# Test of the value obtained
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 #nbrefoissim=
MV[K,3]<-l2/nbrefoissim
MV[K,4]<-l3/nbrefoissim
}
pc
And here is the R code that I used to code mine.
library(MASS)
CentrageV<-function(X,Ms,n){
# this function centres the data from X
X1=X*0
for (i in 1:n){
X1[i,]=X[i,]-Ms
}
return(X1)
}
# Fonction N°2
SqrtMatrice0<-function(M) {
# This function allows us to calculate the square root of a matrix
# using the decomposition M=PDQ where Q is the inverse of P
# here the negative eigenvalues are replaced by 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)
}
# declaration of parameters
m1=0.01 # value of alpha (risque de 1%)
m2=0.05 # value of alpha (risque de 5%)
m3=0.1 # value of alpha (risque de 10%)
nbrefoissim=100 # number of times the program runs
p=3 #dimension of the variable X
q=3 #dimension of the variable Y
R=c(5,4,3);# Number of partitions of each component of the variable Y
if(length(R) != q) stop("The size of R must be equal to q")
n=25 # Taille echantillon
N=c(25,50,100,200,300,400,500,1000) #different sample sizes
#N=c(25,50) #ifferent sample sizes
pc1=rep(0,100)
K=0
MV=matrix(0,nr=8,nc=4)
dimnames(MV)[[2]]=c("n","r1%","r5%","r10%")
#Start of the programme
for (n in N){
l1=0 # initialization of the value to calculate the test level at 1%.
l2=0 # initialization of the value to calculate the test level at 1%. 5%
l3=0 # initialization of the value to calculate the test level at 1%. 10%
# Creation of a list n11 which contains the sizes of the different groups
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])
}
# Creation of the lists P11 and P12 which contain the probabilities and
# the inverses of the empirical probabilities of the different groups
respectively
P11=list()
P12=list()
for (i in 1:q){
P11[[i]]=n11[[i]]/n
P12[[i]]=n/n11[[i]]
}
# creation of a list containing the matrices W
W=list()
for (i in 1:q){
w=matrix(0,n,R[i])
w[1:n11[[i]][1],1]=1
if (R[i]>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){
# data generation
VA1=mvrnorm(n,rep(0,(p+q)),diag((p+q)))
X=VA1[,1:p]
Y=VA1[,(p+1):(p+q)]
# Xbar calculation
Xbar=colMeans(X)
# Calculation of 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
}
#Calculation 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
}
#Calculation 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))
# Calculation of Gamma
GGamma=matrix(0,p*sum(R),p*sum(R))
PGamma=kronecker(diag(P12[[1]]),diag(p))
Ifin=p*R[1]
GGamma[1:Ifin,1:Ifin]=PGamma
for (i in 2:q){
PGamma=kronecker(diag(P12[[i]]),diag(p))
Idebut=((p*sum(R[1:(i-1)]))+1)
Ifin=(p*sum(R[1:i]))
GGamma[Idebut:Ifin,Idebut:Ifin]=PGamma
}
#Calculation of Sigma
# Calculation of Vn
X1=CentrageV(X,Xbar,n)
Vn=t(X1)%*%X1/n
## Calculation of Sigma
GSigma=matrix(0,p*sum(R),p*sum(R))
for (i in 1:q ){
for (j in 1:R[i] ){
Xij=CentrageV(X,Xjhbar[[i]][j,],n)
Vij=t(W[[i]][,j]*Xij)%*%Xij/n11[[i]][j]
for (k in 1:q ){
for (l in 1:R[k]){
Xkl=CentrageV(X,Xjhbar[[k]][l,],n)
Vijkl=t(W[[i]][,j]*Xij)%*%(W[[k]][,l]*Xkl)/n
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))
}
}
}
}
# Eigenvalue determinations of sigma epsilon
pa=SqrtMatrice0(GSigma)
mq= pa %*% GGamma %*% pa
u=Re(eigen(mq)$values)
# determination of degree of freedom and c-value noted va
dl=(sum(u)^2)/(sum(u^2))
va=(sum(u^2))/(sum(u))
pc=1-pchisq(tr1/va, df= dl)
pc1[i1]=pc
# Test of the value obtained
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 #nbrefoissim=number of simulations
MV[K,3]=l2/nbrefoissim
MV[K,4]=l3/nbrefoissim
}
Sincerely
[[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.
