y = X β + ϵ y = X\beta+\epsilon y=Xβ+ϵ y ∈ R n × 1 , X ∈ R n × p , β ∈ R n × 1 , ϵ ∈ R n × 1 y\in \mathbb{R}^{n\times 1},X\in \mathbb{R}^{n\times p},\beta\in \mathbb{R}^{n\times 1},\epsilon \in \mathbb{R}^{n\times 1} y∈Rn×1,X∈Rn×p,β∈Rn×1,ϵ∈Rn×1 数据构造如下 X ∼ N p ( 0 p , Σ ) y = X β + N ( 0 , 0.1 ) X\sim {N_p}({\mathbf{0}_p},{\Sigma })\\ y = X\beta+N(0,0.1) X∼Np(0p,Σ)y=Xβ+N(0,0.1)
生成60个样本，维度为6,协方差为 Σ = d i a g ( 10 , 8 , 6 , 4 , 2 , 1 ) \Sigma=diag(10,8,6,4,2,1) Σ=diag(10,8,6,4,2,1),这样的得到的数据 X X X为满秩，系数可以直接用最小二乘法求解 β ^ = ( X T X ) − 1 X T y \hat{\beta}=(X^TX)^{-1}X^Ty β^=(XTX)−1XTy。
代码主要调用了函数mvnrnd生成多元变量，第一个参数为均值向量，第二个参数为协方差矩阵，第三个参数是生成的样本数量。
n = 60; p = 6; cov = diag([10,8,6,4,2,1]); X = mvnrnd(zeros(1,p),cov,n); %经验协方差估计 hatcov = X'*X/(n-1) beta = rand(p,1); err = mvnrnd(0,0.1,n); y = X*beta + err; hatbeta = inv(X'*X)*X'*y; %画图 plot([beta hatbeta],'linewidth',1.5); h = legend(["$$\beta$$","$$\hat{\beta}$$"]); set(h,'Interpreter','latex') a=gca; a.FontSize=20; a.FontName='Helvetica'; set(gcf,'color','w');协方差估计，与实际的协方差比较接近 系数估计值与实际值基本一致
这里采用了 n ≫ p n\gg p n≫p，如果要生成 n ≪ p n\ll p n≪p要怎么弄呢？
简单的做法如下
X ∼ N p ( 0 p , I p ) X\sim N_p(\mathbf{0}_p,I_p) X∼Np(0p,Ip) 按照这个分布生成 n n n个样本，不过求 β \beta β就不能用上面的公式直接得到，需要用正则或偏最小二乘法等手段，得到稳健一些的 β \beta β
