[R] Help for storing value in the matrix

[Ablaye Ngalaba]

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