In zxy990409/StatComp21015: establish my packages for 'Statistical Computing' course

Overview

StatComp21015 is a simple R package developed to implement two functions which is related to the simulations of the perceived adjacency matrices when we address community detection problems based on the partial network information. Two functions are considered, namely, qua95 (calculating the 0.95 quantile of a new statistic) and subg (generateing the perceived adjacency matrix based on one individual).

The source R code for qua95 is as follows:

qua95 <- function(N){
  rej1 <- list()
  for (z in 1:N) {
    EA <- matrix(0,100,100)
    Q <- matrix(c(0.3,0.1,0.1,0.3), 2, 2)
    z1 <- matrix(c(1,0),2,1)
    z2 <- matrix(c(0,1),2,1)
    EA[1:50, 1:50] <- matrix(t(z1) %*% Q %*% z1, 50, 50)
    EA[51:100, 1:50] <- matrix(t(z2) %*% Q %*% z1, 50, 50)
    EA[1:50, 51:100] <- matrix(t(z1) %*% Q %*% z2, 50, 50)
    EA[51:100, 51:100] <- matrix(t(z2) %*% Q %*% z2, 50, 50)
    A <- matrix(0,100,100)
    for (i in 1:100){
      for(j in 1:100){
        A[i,j] <- rbinom(1,1,EA[i,j])
      }
    }
    A1 <- A
    A1[lower.tri(A1)] <- t(A1)[lower.tri(A1)]
    A3 <- A1[-1,-1]

    b <- A3[1,]
    Nr1 <- length(which (b == 1, arr.ind = T))
    Nr2 <- length(which (b == 0, arr.ind = T))

    u <- (sqrt(Nr1-1)+sqrt(Nr2))^2
    sigma <- sqrt(u)*{(1/sqrt(Nr1-1)+1/sqrt(Nr2))^(1/3)}

    X1 <- matrix(rnorm(Nr1*Nr1,0,1),nrow=Nr1,ncol=Nr1)
    diag(X1) <- rnorm(Nr1,0,2)
    X1[lower.tri(X1)] <- t(X1)[lower.tri(X1)]
    X1 <- X1/sqrt(Nr1)
    X2 <- matrix(rnorm(Nr1*Nr2,0,1),nrow=Nr1,ncol=Nr2)
    rt11 <- abs((max(eigen(t(X2) %*% X2)$values)-u)/sigma)
    rt12 <- abs(max(eigen(X1)$values)-2)*Nr1^(2/3)
    rej1[[z]] <- rt11+rt12
  }

  rej1<- as.vector(unlist(rej1))
  new_rej12 <- rej1[order(rej1)][95]*2 
  return (new_rej12)
}

When we input the number of simulation, we can get the 0.95 quantile of the empirical distribution.

The source R code for subg is as follows:

subg=function(A,L) 
{
  n=dim(A)[1]
  C=B=matrix(rep(0,n^2),nrow=n)
  if (L>1)
  {
    for (i in 1:(n-1))
      for (j in (i+1):n)
        if (A[1,i]==1&A[i,j]==1)
        {B[i,j]=B[j,i]=1
        } else   if (A[1,j]==1&A[i,j]==1)
        {B[i,j]=B[j,i]=1
        } 
  } else {B[1,]=A[1,]
  B[,1]=A[1,]
  }

  if (L>2)
  {
    for (t in 3:L)
    {
      C=matrix(rep(0,n^2),nrow=n)
      for (i in 2:(n-1))
        for (j in (i+1):n)
          if (B[i,j]==1)
          { a=(2:n)[-(i-1)]
          C[i,a]=A[i,a]
          C[a,i]=A[a,i]
          a=(2:n)[-(j-1)]
          C[j,a]=A[j,a]
          C[a,j]=A[a,j]
          }
      B=B+C
      for (i in 1:(n-1))
        for (j in (i+1):n)
          B[i,j]=B[j,i]=min(B[i,j],1)
    }
  }
  return(B)
}

When givein an observed adjacency matrix and the depth of knowledge, we will get the the perceived adjacency matrix, which is useful in the community detection. We want to get 0.95 quantile when the number of sumulation equals to 100, and when we have matrix A and length L=2.

library(StatComp21015)
t <- qua95(100)
A <- matrix(c(0,1,1,1,0,0,1,0,1,0,0,0,1,1,0,0,0,0,1,0,0,0,1,1,0,0,0,1,0,1,0,0,0,1,1,0),6,6)
B <- subg(A,2)
print(t)
print(B)

zxy990409/StatComp21015 documentation built on Dec. 23, 2021, 11:18 p.m.