%% PIII_100Hs（PIII分布）
close all;clear;clc
load test.mat %加载所有数据
lie_number=1; %要提取数据点的列序号
SL=test(:,1); %提取的数据
X=SL'; %矩阵转置
X=sort(X); %sort是把波高从小到大排列
X=fliplr(X); %把波高x从大到小排列
n=length(X); %求出波高资料X序列的长度n
%--计算经验频率，并算出相应的位置坐标，并绘出经验点--
p=[[1:n]/(n+1)]; %求经验频率p
xp=norminv(p,0,1); %求出经验频率p对应的位置坐标xp
x=0:0.001:10; %纵坐标
m=mean(X); %计算样本X的均值
K=X/m; %计算模比系数K
Cv=sqrt((1/(n-1)).*sum((K-1).^2)); %计算变差系数Cv（利用的是《工程水文学》（人民交通出版社 邱大洪主编 第四版）课本51页样本的无偏估计）
Cs=(sum((K-1).^3))/(Cv^3)/(n-3); %计算偏态系数Cs
alfa=4/(Cs^2); %P3分布参数计算
beta=(m*Cv*Cs)/2; %P3分布参数计算
a0=m-2*Cv*m/Cs; %P3分布参数计算
P=1-gamcdf(x-a0,alfa,beta); %计算理论频率P
XP=norminv(P,0,1); %计算理论频率P所对应的坐标XP
% p3分布 CDF图--------------------------------------------------------------
% cumulative distribution function plots empirical and modeled 经验累积分布图cdf%
figure1 = figure; %准备画图
axes1 = axes('Parent',figure1); %坐标轴
plot(sort(SL,'ascend'),p,'linewidth',2); %实际数据
hold on %在原图上继续添加内容
plot(sort(x,'ascend'),1-P,'r','linewidth',2) %用P3拟合的数据
title('CDF plot') %设定标题
set(gca,'FontWeight','bold'); %设定字体
xlabel('Wave Height(m)','FontName','Times New Roman','FontWeight','bold'); %设定x轴标题
ylabel('Cumulative Frequency','FontName','Times New Roman','FontWeight','bold'); %设定y轴标题
backColor = [0.68627 0.93333 0.93333]; %设定背景色的RGB值
set(gca, 'color', backColor) %设定背景色
% print(gcf,'-dtiff',strcat(['P3分布CDF图第',num2str(lie_number),'列.tif']));%保存图片
% ---------------------------------------------------------------------------------
% p3分布 Q-Q图--------------------------------------------------------------
figure(2)
gqqplot(SL,'gam') %直接用evim工具包中的函数画图
title('Q-Q plot') %设定标题
set(gca,'FontWeight','bold'); %设定字体
xlabel('Theoretical Quantiles','FontName','Times New Roman','FontWeight','bold'); %设定x轴标题
ylabel('Quantiles of Input Sample','FontName','Times New Roman','FontWeight','bold'); %设定y轴标题
backColor = [0.68627 0.93333 0.93333]; %设定背景色的RGB值
set(gca, 'color', backColor) %设定背景色
% print(gcf,'-dtiff',strcat(['P3分布Q-Q图第',num2str(lie_number),'列.tif']));%保存图片
% ---------------------------------------------------------------------------------
figure(3)
%横坐标：x=norminv(p,mu,sigma),标准正态分布时，mu=0,sigma=1)
plot(xp,X,'.','markersize',8) %绘出经验点
hold on %在原图上继续添加内容
%--计算经验频率，并算出相应的位置坐标，并绘出经验点--
plot(XP,x,'r-.','LineWidth',1.6) %绘出拟合曲线
PL=[0.01 0.05 0.5 1 2 5 10 20 30 40 50 60 70 80 90 95 98 99 99.9 99.99];%设置海森机率格纸的坐标轴,PL为频率
xx=norminv(PL/100,0,1); %PL所对应的海森机率格纸的横坐标
yy=[0:0.5:5]; %纵坐标的间隔设置
axis([min(xx),max(xx),floor(X(end)),ceil(X(2))]) %设置坐标轴的范围
set(gca,'xtick',[xx],'ytick',[yy]) %设置坐标轴的刻度
set(gca,'xticklabel',[PL],'yticklabel',[yy]) %设置坐标轴的刻度
grid on %绘制海森机率格纸
title('P3 distribution'); %设定标题
xlabel('Frequency(%)','FontName','Times New Roman','FontWeight','bold');%设定x轴标题
ylabel('H(m)','FontName','Times New Roman','FontWeight','bold'); %设定y轴标题
% print(gcf,'-dtiff',strcat(['P3分布图第',num2str(lie_number),'列.tif']));%保存图片
%计算不同重现期的重现波高，更改pper的值，100年即pper=100，50年即pper=50--------
%Return100为例
pper=100; %要计算的重现期
%利用try函数找出100年重现期对应的累积频率0.01，误差为0.0001
try
allP=(find(P>1/pper-0.0001 & P<1/pper+0.0001));
Hreturn=x(allP(floor(length(allP)/2)))
catch
allP=(find(P>1/pper-0.001 & P<1/pper+0.001));
Hreturn=x(allP(floor(length(allP)/2)))
end
%Hreturn即为对应重现期的重现波高
评论1