基于MASS包的mvrnorm产生多元高斯分布的随机数

太阳不下山

library(MASS)
library(Matrix)

定义一些常数

b <- 10
q <- 10
p <- 100
n <- 50
r <- 0.95

根据说明，X（nXp）是产生自10个独立同分布的多元高斯随机分布(10组随机变量，每组随机变量高度相关) 即每组随机变量服从N（0，Sigma), Sigma(i,i)=1, Sigma(j,k)=0.95 (j!=k) 定义函数，生成X

generate_X <- function(q,b,n,r){
  # q: 每组随机变量个数
  # b: 同分布的随机变量组数
  # n: 随机样本大小
  # r: 同组随机变量间的任意两个随机变量的相关系数
 
   mu <- rep(0,q*b) # 多元高斯分布的均值向量（均值都为0）# 100维
    Sigma <- matrix(r, ncol=q,nrow=q) # 初始化协方差阵为全部元素0.95的矩阵
    diag(Sigma) <- 1 # 协方差阵的对角线更正为1

基于MASS包的mvrnorm产生多元高斯分布的随机数，由于每组随机变量同分布， 故产生b组随机数，即为X一行，重复n次，即得到随机样本X

Sigmas <- as.matrix(bdiag(lapply(1:b,function(o) Sigma))) # 生成100维X100维大协方差矩阵，由10个Sigma组成的对角矩阵块
  X <- mvrnorm(n=n,mu=mu, Sigma=Sigmas) # 产生服从N（0，Sigmas)的随机数
  return(X)
}

定义函数，生成y

generate_y <- function(X){
  # X:设计矩阵
  # y = -X1+2X11-3X21+4X31-5X41+6X51-7X61+8X71-9X81+10X91+E
  # 根据上面公式提取出beta系数
  beta <- rep(0,ncol(X)) # 初始化beta系数为与X相同列数的0向量
  idx <- seq(1,100,10) # 生成1,11,21,31,...,91
  beta[idx] <- c(-1,2,-3,4,-5,6,-7,8,-9,10) # 更正每组第一个变量的系数
  y <- X%*%beta+rnorm(nrow(X)) # X.beta+E 即为得到y，E为标准正态随机数（rnorm()生成）
  return(y)
}

分别生成训练集、测试集、验证集 training set

#train_x <- generate_X(q,b,n,r)
#train_y <- generate_y(train_x)

testing set

#test_x <- generate_X(q,b,n,r)
#test_y <- generate_y(test_x)

validation set

#valid_x <- generate_X(q,b,n,r)
#valid_y <- generate_y(valid_x)
 x=generate_X(q,b,n,r);
  y=generate_y(x);
 p=1/(1+exp(-y))
 Y= rep(0,nrow(p)); for(i in 1:nrow(p)){Y[i]=rbinom(1,1,p[i])}
 Y= as.matrix(Y)
 data=cbind(Y,x)