p_iii_new_3.m

皮尔逊三型算法的源程序,比较好用.希望大家都能用得好.

% "P-III"
clear all;
close all
%fid=fopen('gumb.dat','r');
%R='pearson-iii.txt';
%R='zhuhai.txt';
R='e:\tide.txt';
fid=fopen(R,'r');

Sta=fscanf(fid,'%s\n',1);
nn=fscanf(fid,'%d\n',1);
for i=1:nn
%fprintf(fid,'%10s\n',leg(i).l);
sta(i).s=fscanf(fid,'%s',1);
zm(i).s=fscanf(fid,'%s\n',1);
mm(i)=fscanf(fid,'%d\n',1);
Y(i).y=fscanf(fid,'%f\n',mm(i));
%fscanf(fid,'\n');
end
%fid=fopen('wave.dat','r');
%fid=fopen('wv.dat','r');
%[H,count]=fscanf(fid,'%f',inf);
fclose(fid);
display('No. of Stations  =');
disp(nn);
n0=1;
%for n0=1:4
while n0~=999
n0=input('Select No. of Station:(if input 999 then STOP)');
if n0==999
    break
end
H=Y(n0).y;
%Hs = load('atlantic.dat');
n=mm(n0);
H=flipud(sort(H));
Hs=H;
Hm=mean(Hs);
ki=Hs/Hm;
Cv1=std(H)/Hm;
Cv=sqrt(sum((ki-1).^2)/(n-1));Cv1=std(Hs)/Hm;
Cs=sum((ki-1).^3)/(n-3)/Cv^3;
%抽样误差各统计参数均方误公式
%P=empdistr(H);
P=[1:n]/(n+1);
PP=[.01,0.1,0.2,0.5,1,2,5,10,20,50,75,90,95,99];
m=length(PP);
display(Cv);display(Cs);
Hs=flipud(Hs);
coef(1)=1;coef(2)=2;coef(3)=3;coef(4)=4;coef(5)=5;coef(6)=6;
A0=zm(n0).s
A1='  测站:';
STA=[A1,A0];
for k=1:6
    if k~=1
       Cs=(coef(k))*Cv;
   end       
a=4/Cs^2;
b=2/(Hm*Cv*Cs);
a0=Hm*(1-2*Cv/Cs);
Tp=a0+gaminv(1-PP/100,a,1)/b;
%Hp0=Tp/b+a0;
Hp0=xp(PP/100,Cs,b,a0);
Hp(:,k)=Hp0';
Tr=1./PP*100;

format short e
data=[H,P'];varin='.';
Cv0=num2str(Cv,3);Cs0=num2str(Cs,3);
leg1='Cv=';leg2=num2str(Cv0);leg3='Cs=';leg4=num2str(coef(k));leg5='Cv';if k==1;leg4=Cs0;leg5='';end
%leg5='\timesCv';
leg=[leg1,leg2,'  ',leg3,leg4,leg5];
ex1=char('  经验累积频率分布');
ex=char('  理论累积频率曲线',STA,leg);
figure(k)
set(gcf,'Color','w');
%h_ap=axes('Position',[0.13,0.13,0.8,0.75]);
normprob_plot_1(data,varin,Hp(:,k),PP/100);
legend(ex1,ex);
%set(gca,'FontSize',10,'FontWeight','bold','fontname','宋体')
%xlabel('水位 (m^3/s)','FontSize',12,'FontWeight','bold','fontname','宋体');
%set(h_ap,'Yaxislocation','left')						%<15>
%ylabel('降雨量 (mm)','FontSize',12,'FontWeight','bold','fontname','宋体');
%set(h_ap,'Xaxislocation','bottom')										%<14>
%xlabel('P(%)','FontSize',12,'FontWeight','bold','fontname','new norman');
%title('P-III型分布','fontname','宋体')
%text(1.6,log(0.03),leg,'FontSize',10,'FontWeight','demi','FontName','宋体');
%grid on
%set(gca,'GridLineStyle','-');


P0= interp1(Hp(:,k),PP'/100,H,'spline'); %specifies alternate methods.

Q(k)=sum(P0-P').^2;
  HH=P0-P';
alpha=0.05;
[h,p,ks] = kstest2(P0,P',alpha,0);
[H0,pval]=chi2test(HH', alpha);

%fid = fopen('d:\shuiwen\P-III-2.txt','a+');
fprintf(1,'%30s\n',R);
fprintf(1,'%12s \n',STA);
fprintf(1,'%12s  %12s\n','Cv','Cs');
fprintf(1,'%12.3f  %12.3f\n',Cv,Cs);
Out=[PP',Tr',Hp(:,k)];
fprintf(1,'%12s  %12s  %12s\n','频率/%','重现期/yr','计算值');
fprintf(1,'%12.4f  %12.2f  %12.2f\n',Out');
fprintf(1,'%10s  %10s  %10s  %10s  %10s  %10s  %10s\n','H-KS假设','P-Value','KS统计值','Alpha','H-卡方假设','P-value','Q统计值');
fprintf(1,'%10i  %10.3f  %10.3f  %10.2f  %10i  %10.3f  %10.3f\n',h,p,ks,alpha,H0,pval,Q(k));
%fclose (fid)

%nam1='d:\shuiwen\fig\Piii';nam2=num2str(k);nam3='.jpg';
nam1='d:\shuiwen\fig\Piii';nam2=num2str(k);nam3='.png';
%nam1='d:\shuiwen\fig\Piii';nam2=num2str(k);nam3='.tif';
nam=[nam1,nam2,nam3];
[X,map]=getframe(k);
imwrite(X,map,nam);


%k=k+1;
end
%clear H;clear data;clear Hp;clear PP;
end