超立方采样程序(Cholesky分解)

基于拉丁超立方采样的MonteCarlo概率潮流计算方法对具有相关性的输入随机变量进行超立方采样的过程，最终生成采样矩阵S。Matlab代码：clc;clearall;Cx1=[1.00.50.80.70.6;0.51.00.70.60.8;0.80.71.00.50.6;0.70.60.51.00.9;0.60.80.60.91.0];%正定矩阵的Cholesky分解[m,n]=size(Cx1);%[5,5]ifm~=n%判断输入的矩阵是不是方阵disp('输入的矩阵不是方阵，请重新输入');return;endfori=1:n%判断输入的矩阵是不是对称矩阵forj=1:nifCx1(i,j)~=Cx1(j,i)disp('输入的方阵不是对称矩阵，请重新输入');return;endendendd=eig(Cx1);%根据方阵的特征值判定是不是正定矩阵fori=1:nifd(i)==0disp('输入的矩阵不是正定矩阵，请重新输入');return;elsebreak;endenddisp('输入的矩阵可以进行Cholesky分解');%如果是正定矩阵，可以进行下面的分解操作%Cholesky分解下三角矩阵R1=chol(Cx1,'lower');randn('state',sum(100*clock));%利用时钟设置随机种子，这样每次产生的随机数就不同了N=1000;%1.设置样本个数W1=zeros(5,N);WW1=zeros(5,N);fori=1:NW1(:,i)=randn(5,1);%5行1列的随机数%W1是符合标准正态分布变量样本endWW1(1,:)=W1(1,:).*0.45+4.5;%产生一个高斯分布的指定均值和方差的矩阵：将randn产生的结果乘以标准差，然后加上期望均值即可WW1(2,:)=W1(2,:).*0.63+3.7;WW1(3,:)=W1(3,:).*0.58+5.3;WW1(4,:)=W1(4,:).*0.51+4.9;WW1(5,:)=W1(5,:).*0.56+4.2;Z1=R1*W1;rouMW1=corrcoef(W1');%转置矩阵的相关系数rouMZ1=corrcoef(Z1');%2.生成顺序矩阵：对矩阵中每行中的某数在该行中大小排名值，即在此行数值大小的中排第几即为几Ls1=zeros(5,N);forp=1:5maxZ1=max(Z1(p,:));%每行最大值赋值maxZ1k=1;forq=1:N[minZlocZ]=min(Z1(p,:));%输出每行中最小值及列数Ls1(p,locZ)=k;k=k+1;Z1(p,locZ)=maxZ1+1;endendLLs1=Ls1;%3.需要进行抽样的的实际变量x1=zeros(5,N);fori=1:Nsx1=(i-0.5)/N;x1(1,i)=norminv(sx1,4.5,0.45);%正态分布μ=4.5，σ=0.45x1(2,i)=norminv(sx1,3.7,0.63);x1(3,i)=norminv(sx1,5.3,0.58);x1(4,i)=norminv(sx1,4.9,0.51);x1(5,i)=norminv(sx1,4.2,0.56);end%4.根据顺序矩阵生成超立方抽样后的样本S1=zeros(5,N);fori=1:5tx1=x1(i,:);tLs1=LLs1(i,:);maxy1=max(tx1);maxtLs1=max(tLs1);forp=1:N[C1D1]=min(tx1);[C2D2]=min(LLs1(i,:));tLs1(1,D2)=C1;tx1(1,D1)=maxy1+1;LLs1(i,D2)=maxtLs1+1;endS1(i,:)=tLs1;endrouMS1=corrcoef(S1');