基于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 <- 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=rbinom(1,1,p)}
Y= as.matrix(Y)
data=cbind(Y,x)
页:
[1]