Mix-Copula.zip_Mix-Copula_copula 注释_mix_混合Copula函数估计_辅助学习

clear clc % N=1000; alpha1=0.6; alpha2=0.4; N1=alpha1*N; N2=alpha2*N; rho1=0.3; rho2=0.9; xx1=copularnd('Gaussian',rho1,N1); xx2=copularnd('Gaussian',rho2,N2); xx=[xx1;xx2]; xxx=xx(:,1); yyy=xx(:,2); U=xxx; V=yyy; %利用EM算法进行估计运算 S=2; % S为Copula个数 N=length(U); % N为样本容量 th=0.000001; %收敛条件 a=2.21; gama=0.65; %参数初始化 rana=[1/2,1/2]'; %权重参数初始化 a1=copulafit('Gaussian',[U V]); a2=copulafit('Gaussian',[U,V])+0.1; theta=[a1(1,2),a2(1,2)]'; %相依结构参数初始化 t=inf; count=0; %迭代次数 COPULA=zeros(N,S); %COPULA函数矩阵 ranack=zeros(S,1); %λ(k)*copula(n,k) n固定 RZNK=zeros(N,S); %λ(k)*copula(n,k)/sum（λ(k)*copula(n,k)） SCADD=zeros(S,1); %λ(k)*SCAD(λ(k)) turntheta_old=zeros(S,1); linss1=zeros(2,51); %用来存储各个Copula权重参数的变化 linss2=zeros(2,51); %用来存储各个Copula相依参数的变化 while t>=th&&count<=50 rana_old=rana; theta_old=theta; turntheta_old(1)=tan(pi*theta_old(1)/2); turntheta_old(2)=tan(pi*theta_old(2)/2); %由于相依结构参数有约束限制，故将参数进行转换，从而用无约束最优化方法来解决 %计算COPULA矩阵 for n=1:N COPULA(n,1)=copulapdf('Gaussian',[U(n,:) V(n,:)],theta_old(1)); %Gaussian Copula密度函数值 COPULA(n,2)=copulapdf('Gaussian',[U(n,:) V(n,:)],theta_old(2)); %Gaussian Copula密度函数值 end %计算SCADD for k=1:S SCADD(k)=rana_old(k)*SCAD(rana_old(k),a,gama); end delta=N*(1-sum(SCADD)); %E步 %开始求rznk for cn=1:N for cm=1:S ranack(cm)=rana_old(cm)*COPULA(cn,cm); end for cm=1:S RZNK(cn,cm)=ranack(cm)/sum(ranack); end end %M步 %求权重参数的更新方程 for cm=1:S rana(cm)=(sum(RZNK(:,cm))-N*rana_old(cm)*SCAD(rana_old(cm),a,gama))/delta; end if sum(rana<0)>=1 rana=rana_old; break;%判断如果出现rana小于0的情况，则终止跳出循环 end rana; linss1(:,count+1)=rana; %求相依结构参数的更新方程 %采用BFGS方法求相依结构参数的更新方程 [turntheta,fval,exitflag,output,grad,hessian]=runtheta(N,U,V,RZNK,t); turntheta; %再把相依结构参数转换回来 theta(1)=atan(turntheta(1))*2/pi; theta(2)=atan(turntheta(2))*2/pi; theta; linss2(:,count+1)=theta; %判断终止条件 t=max([norm(rana_old(:)-rana(:));norm(theta_old(:)-theta(:))]); count=count+1; end count = 51 theta theta = 0.1998 0.8889 rana rana = 0.5920 0.4080 theta = 0.2419 0.8898 rana rana = 0.5404 0.4596