利用least-square参数估计法与Copula函数计算联合概率分布 (matlab code)

%利用least-square参数估计法与Copula函数计算联合概率分布
%parameter estimation for Weibull distribution and estimate the parameter for Copula
clear;clc;
n=20;
h=[8.7 12.8 9.9 10.3 7.5 12.3 10.3 7.5 9.3 12.7 6.8 7.4 12.6 12.2 8.7 13.2 11.8 11.0 10.5 4.6];
s=[39.9 40.3 41.7 28.5 26.4 47.1 40.5 22.9 40.5 42.2 23.8 23.6 38.2 45.8 34.2 31.2 35.2 41.2 41.5 22.4];
%--------determine Weibull distribution-----------------------------------
h_sort=sort(h);
s_sort=sort(s);
x_hi=log(h_sort);
x_si=log(s_sort);
for k=1:20
MR(k)=(k-0.3)/(n+0.4);
x_hi_2(k)=(x_hi(k))^2;
x_si_2(k)=(x_si(k))^2;
end
yi=log(-log(1-MR));
for i=1:n
yi_2(i)=(yi(i))^2;
xy_hi(i)=x_hi(i)*yi(i);
xy_si(i)=x_si(i)*yi(i);
end
b_hat_h=(sum(xy_hi)-sum(x_hi)*sum(yi)/n)/(sum(x_hi_2)-(sum(x_hi))^2/n);
b_hat_s=(sum(xy_si)-sum(x_si)*sum(yi)/n)/(sum(x_si_2)-(sum(x_si))^2/n);
a_hat_h=mean(yi)-b_hat_h*mean(x_hi);
a_hat_s=mean(yi)-b_hat_s*mean(x_si);
eta_hat_h=exp(-a_hat_h/b_hat_h)
eta_hat_s=exp(-a_hat_s/b_hat_s)
beta_hat_h=b_hat_h
beta_hat_s=b_hat_s
%----------determine Copula function--------------------------------------------
count_nag=0;
count_pos=0;
for j=1:n-1
for l=j+1:n
tau_value=(h(j)-h(l))*(s(j)-s(l));
if(tau_value<0)
count_nag=count_nag+1
else
count_pos=count_pos+1;
end
end
end
tau=(count_pos-count_nag)/(n*(n-1)*0.5)
theta=2*tau/(1-tau)

%Plot pdf and cdf distribution figures in calculation
%Prepare by Pu Li, 2008.03.22
clear;clc;
h=[0:0.3:30];
s=[0:0.9:90];
h_pdf=wblpdf(h, 11.0175, 4.2393);
s_pdf=wblpdf(s, 38.8137, 4.4892);
figure;
plot(h, h_pdf);
figure;
plot(s, s_pdf);
h_cdf=wblcdf(h, 11.0175, 4.2393);
s_cdf=wblcdf(s, 38.8137, 4.4892);
figure;
plot(h, h_cdf);
figure;
plot(s, s_cdf);
for i=1:101
for j=1:101
Copula_C(j)=(((h_pdf(i)))^(-1.9583)+((s_pdf(j)))^(-1.9583)-1)^(-1/1.9583);
Copula_G(j)=exp(-((-(log(h_pdf(i))))^1.9583+(-(log(s_pdf(j))))^1.9583)^(1/1.9583));
Copula_CC(j)=(((h_cdf(i)))^(-1.9583)+((s_cdf(j)))^(-1.9583)-1)^(-1/1.9583);
Copula_GC(j)=exp(-((-(log(h_cdf(i))))^1.9583+(-(log(s_cdf(j))))^1.9583)^(1/1.9583));
end
sub_cc(i,:)=real(Copula_C);
sub_cg(i,:)=real(Copula_G);
sub_ccc(i,:)=real(Copula_CC);
sub_cgc(i,:)=real(Copula_GC);
end
figure;
surf(s, h, sub_cc);
figure;
surf(s, h, sub_cg);
figure;
surf(s, h, sub_ccc);
figure;
surf(s, h, sub_cgc);

%Model fitting test
%Prepare by Pu Li, 2008.03.22
clear;clc;
n=20;
h=[8.7 12.8 9.9 10.3 7.5 12.3 10.3 7.5 9.3 12.7 6.8 7.4 12.6 12.2 8.7 13.2 11.8 11.0 10.5 4.6];
s=[39.9 40.3 41.7 28.5 26.4 47.1 40.5 22.9 40.5 42.2 23.8 23.6 38.2 45.8 34.2 31.2 35.2 41.2 41.5 22.4];
%-----------------Calculate the experimental cdf-------------------------------------------
for i=1:n
n_sample=0;
for j=1:n
if(h(j)<=h(i) & s(j)<=s(i))
n_sample=n_sample+1;

end
n_sample_collect(i)=n_sample;

end
%-------------------------Calculate Theoritical cdf-----------------------------
cp_h(i)=1-exp(-(h(i)/11.0175)^4.2393);
cp_s(i)=1-exp(-(s(i)/38.8137)^4.4892);
Copula_C(i)=(((cp_h(i)))^(-1.9583)+((cp_s(i)))^(-1.9583)-1)^(-1/1.9583);
Copula_G(i)=exp(-((-(log(cp_h(i))))^1.9583+(-(log(cp_s(i))))^1.9583)^(1/1.9583));
end
%----------------------------------------------------------------------
P_sample=n_sample_collect/n;
sub_cc=real(Copula_C);
sub_cg=real(Copula_G);
%----plot the theoritical vs experimental figures and calculate the correlation coefficient--
figure;
plot(sub_cc, P_sample,'*');
figure;
plot(sub_cg, P_sample,'*');
Linear_coef_cc=corrcoef(sub_cc, P_sample)
Linear_coef_cg=corrcoef(sub_cg, P_sample)