基于单变量和多变量计算水文频率-基础知识&程序实现

内容一：概率论基础
随机变量分布函数
举例：某水文站的水位
分布函数：设X是一个随机变量，对任意实数x，称F(x)=P（X≤x）为随机变量X的分布函数。
水文频率：人们根据水文现象的随机性，用概率来描述未来出现各种大小洪水的可能性。设计时，对给定的概率p（水文上称为频率），选择满足关系式1-1的Xp 作为工程的设计依据。
P（X>=Xp）=p
p称为设计频率，Xp称为设计值。重现期T=1/P。
例如：某河段的洪峰流量超过3000m3/s时的概率为1%，则设计值为3000m3/s, 依据此建造的河段闸坝设计频率为1%，设计标准为100年一遇。
水文频率计算的目的：给定设计频率P，求对应设计值Xp。
内容二：基于单变量计算水文频率
分布线型：
PIII，GEV，LN2
参数估计：
极大似然法，L矩法，
拟合优度检验
KS检验，AIC，BIC准则，Rmse误差
附：Matlab代码
内容三：基于多变量计算水文频率
分布线型
参数估计
拟合优度检验
附：Matlab代码

程序实现：（新手上路，代码有不足的地方还请路过的高手批评指正）
lear all;close all;clc;
data=xlsread(‘copula2.xlsx’);%%读取原始数据
%% 参数估计
for i=1:5;
idata=data(:,i);%读取站点数据列
% 读取P3分布参数
parap3=xlsread(‘para.xlsx’);
iEx=parap3(i,1);
iCv=parap3(i,2);
iCs=parap3(i,3);
%GEV分布参数
[iparmhatgev,iparmcigev]=gevfit(idata);
ik=iparmhatgev(1);
isigma=iparmhatgev(2);
imu=iparmhatgev(3);
% LN2分布参数
[iparmhatlogn,parmcilogn]=lognfit(idata);
ilnavg=iparmhatlogn(1);
ilnsd=iparmhatlogn(2);
% 将3种分布的参数显示到表格Result中
Result(i,1)=iEx;
Result(i,2)=iCv;
Result(i,3)=iCs;
Result(i,4)=ik;
Result(i,5)=isigma;
Result(i,6)=imu;
Result(i,7)=ilnavg;
Result(i,8)=ilnsd;
%% KS检验之P3
n=length(idata);
idatasort=sort(idata);%从小到大排列
real_data=sort(idata,‘descend’);
for m=1:n;
iexper(m)=length(find(idatasort<=idatasort(m)))/n;%实测经验分布，累积分布函数
iex(m)=length(find(real_data>=real_data(m)))/(n+1);
% syms x %定义ix为符号变量
ialpha=4/iCs^2;
% ibeta=2/(iExiCviCs);
% ia0=iEx*(1-2iCv/iCs);
% iy1=((x-ia0)^(ialpha-1))(exp(-ibeta*(x-ia0)));
% ify1(m)=int(iy1,x,idatasort(m),+inf);
% ip(m)=(ibeta^ialpha)/gamma(ialpha)ify1(m);
% ipp(m)=1-ip(m);%%pp是PIII分布的累积分布函数
end
% ippp=double(ipp);
iexperience=iexper’;
% itheory=ippp’;
iexx=iex’;
imean=mean(idata);
itheory2=1-gamcdf((2(real_data-imean)/(imeaniCviCs)+ialpha),ialpha,1);
[iasiran,ikssiran]=kstest(iexx,[iexx,itheory2],0.05);%%KS检验 yy是经验频率% a_siran=1时，拒绝原假设，=0时，接受原假设
ikssi=max(abs(itheory2-iexx));%KS检验D值 %abs去绝对值
iAICsi=nlog(sum((itheory2-iexx).^2)/n)+22;%AIC准则
iOSLsi=sqrt(sum((itheory2-iexx).^2)/n);%OSL准则，即离差平方和最小准则
%% KS检验之GEV分布
igevcdf=gevcdf(idatasort,ik,isigma,imu);%%累积分布函数,理论频率，由大到小（由小值概率到大值概率）
%iyy=(1:length(DS))’./(length(DS)+1);%此经验频率算出来是由小到大的（由大值的概率到小值的概率）
%iyyy=sort(iyy,‘descend’);%经验频率，重新排列后成为由小值概率到大值概率
[ilogicalgev,icriticalgev]=kstest(iexperience,[iexperience,igevcdf],0.05);%%KS检验 yy是经验频率% a_siran=1时，拒绝原假设，=0时，接受原假设
iDgev=max(abs(igevcdf-iexperience));%KS检验D值 %abs去绝对值
iAICgev=nlog(sum((igevcdf-iexperience).^2)/n)+22;%AIC准则
iOSLgev=sqrt(sum((igevcdf-iexperience).^2)/n);%OSL准则，即离差平方和最小准则
%% KS检验之Lognormal分布
ilogncdf=logncdf(idatasort,ilnavg,ilnsd);%%累积分布函数,理论频率，由大到小（由小值概率到大值概率）
%iyy=(1:length(DS))’./(length(DS)+1);%此经验频率算出来是由小到大的（由大值的概率到小值的概率）
%iyyy1=sort(iyy,‘descend’);%经验频率，重新排列后成为由小值概率到大值概率
[ilogicallgn,icriticallgn]=kstest(iexperience,[iexperience,ilogncdf],0.05);%%KS检验 yy是经验频率% a_siran=1时，拒绝原假设，=0时，接受原假设
iDlgn=max(abs(ilogncdf-iexperience));%KS检验D值 %abs去绝对值
iAIClgn=nlog(sum((ilogncdf-iexperience).^2)/n)+22;%AIC准则
iOSLlgn=sqrt(sum((ilogncdf-iexperience).^2)/n);%OSL准则，即离差平方和最小准则
%% KS检验结果体现在表格Result2中
Result3(i,1)=iasiran;
Result3(i,2)=ilogicalgev;
Result3(i,3)=ilogicallgn;
Result2(i,2)=ikssiran;
Result2(i,3)=ikssi;
Result2(i,4)=iAICsi;
Result2(i,5)=iOSLsi;
Result2(i,7)=icriticalgev;
Result2(i,8)=iDgev;
Result2(i,9)=iAICgev;
Result2(i,10)=iOSLgev;
Result2(i,12)=icriticallgn;
Result2(i,13)=iDlgn;
Result2(i,14)=iAIClgn;
Result2(i,15)=iOSLlgn;
end