利用iml模拟多元正态分布

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

前阵子，做了一些纵向数据的模拟，其中都牵涉到
多元正态分布的模拟。
与一元分布相比，多元分布主要加入了不同随机变量
之间的相关关系，即协方差矩阵。
这里介绍一下，以SAS/iml生成多元正态分布数据；
iml中提供直接生成多元正态分布的函数，
当然我们也可以根据多元正态的定义，来生成随机序列，
其结果是一样的。
proc iml;
call randseed(1);/*设定随机数的seed,方便下面比较*/
N=5;             /*随机序列的数目，即观测数*/
mu={1 2 3};      /*多元正态分布的均值*/
cov={4 1 2,
     1 9 -3,
     2 -3 25};
res=randnormal(n,mu,cov);/*通过randnormal函数，生成多元正态分布序列*/

/*作为比较，下面通过多元正态分布的定义，来生成序列*/

std_n=j(N,3,0);
call randseed(1,1);  /*令seed重置，即重复上述的随机数序列*/
call randgen(std_n,'normal',0,1); /*将std_n变为标准正态分布序列，
                                    此时，列与列之间是相互独立的*/
sig=root(cov);                    /*对协方差矩阵进行cholesky分解，生成非奇异
                                    上三角矩阵矩阵sig,且sig'*sig=cov */
res1=std_n*sig+repeat(mu,N,1);    /*通过定义生成多元正态分布序列*/
print res res1;                   /*res与res1是一样的*/
quit;

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏6 回帖

关键词：多元正态分布 正态分布 Cholesky normal choles 正态分布

/* my method for multivariate rand normal dist*/

data a (type=corr);
_type_='corr';
input x1-x3;
cards;
1   .   .
.3  1   .
.5 .4   1
;
proc factor n=3;
run;

proc iml;
f={    0.78114        -0.46758         0.41375,
    0.70560         0.68785         0.17028,
    0.83518        -0.14380        -0.53084 };
data=rannor(j(5000,3,0));
z = data * t(f);     /*impose inner correlation*/
x1=z[,1]*.1 + 5; /*linear transformation*/
x2=z[,2]*.2 + 3;
x3=z[,3];
z=x1||x2||x3;
create dta1 from z [colname={x1 x2 x3}];
append from z;
quit;

proc corr data=dta1;
var x1 x2 x3;
run;

davil2000 发表于 2012-12-17 15:03
/* my method for multivariate rand normal dist*/

data a (type=corr);

首先谢谢davil2000；
这个方法我在网上也看到过；
通过相关阵，利用因子分析的过程求得因子荷载；
这里
proc factor 后面没有其他选项，即用默认的主成分法；
主成分法，从矩阵角度来看，
即对相关阵进行了特征值分解，
假设这里的特征值为e1,e2,e3（e1>=e2>=e3）;
对应的特征向量为t1,t2,t3（均为列向量）;
而后的因子荷载为[sqrt(e1)*t1,sqrt(e2)*t2,sqrt(e3)*t3];
按照因子的分析的算法：
proc iml;
a={1   .3   .5,
  .3  1   .4,
  .5 .4   1};
b=eigval(a);
c=eigvec(a);
print b c;
d=repeat(sqrt(b`),3,1)#c;
print d;
quit;
求得的d也与proc factor的因子荷载相等。
不过我不理解的是：
data=rannor(j(5000,3,0));
z = data * t(f);     /*impose inner correlation*/
这一步来给变量之间加相关关系，总觉得与多元正态的
定义不太符合，求解释~

此外，proc factor，
好像默认会使用相关阵进行计算，即使输入协方差阵，也会先
转为相关阵再计算。
从协方差到相关阵，还是丢掉了变量本身的方差信息。

ziyenano 发表于 2012-12-17 16:15
首先谢谢davil2000；
这个方法我在网上也看到过；
通过相关阵，利用因子分析的过程求得因子荷载；

“总觉得与多元正态的定义不太符合…”？
正态分布的变量(们)经线性组合后产生的新的变量仍然是正态分布。这种新变量当然与组成的元素变量之间存在相关关系了。
至于协方差阵转换到相关阵过程中会产生信息损失的说法，解释为统计指标都会有信息损失，无信息损失的只能是原始数据了。O(∩_∩)O~

ziyenano 发表于 2012-12-17 16:15
首先谢谢davil2000；
这个方法我在网上也看到过；
通过相关阵，利用因子分析的过程求得因子荷载；

“总觉得与多元正态的定义不太符合…”？
正态分布的变量(们)经线性组合后产生的新的变量仍然是正态分布。这种新变量当然与组成的元素变量之间存在相关关系了。
至于协方差阵转换到相关阵过程中会产生信息损失的说法，解释为统计指标都会有信息损失，无信息损失的只能是原始数据了。O(∩_∩)O~

davil2000 发表于 2012-12-17 16:31
“总觉得与多元正态的定义不太符合…”？
正态分布的变量(们)经线性组合后产生的新的变量仍然是正态分布 ...

怪我问题问得不好~
其实我的疑问是，如何保证以因子荷载的矩阵，对变量进行线性变换，
使得变量之间的相关性恰好为开始设定的相关阵；
这种方法的结果，我以前测验过，是没有问题的。
莫非以[sqrt(e1)*t1,sqrt(e2)*t2,sqrt(e3)*t3]，对变量进行线性变化，
得到的变量之间的相关性，一定是服从原始相关阵的？
我觉得，数学上应该可以推导的。

不过要想以协方差阵构造多元正态，这个方法就有点烦了。

很有启示，已收藏！