import numpy as np from math import * import sys from scipy.stats import gamma,norm from scipy.special import gamma as Gamma import matplotlib.pyplot as plt def NEP(n,alpha): ii=np.arange(1,n+1,1) return np.array((ii-alpha)/(float(n+1)-2.0*alpha)) def MSR(n,z): # mz=np.average(z) # sz=np.std(z,ddof=1) # _s=np.std(z,ddof=0) mz=np.sum(z)/float(n) sz=np.sqrt(np.sum((z-mz)**2)/float(n-1)) _s=np.sqrt(np.sum((z-mz)**2)/float(n)) _c=np.sum(((z-mz)/_s)**3)/float(n) rz=sqrt(float(n)*float(n-1))/float(n-2)*_c return mz,sz,rz def LMOM(n,z): jj=np.arange(1,n+1,1) b0=np.sum(z)/float(n) b1=np.sum((jj-1)*z)/float(n)/float(n-1) b2=np.sum((jj-1)*(jj-2)*z)/float(n)/float(n-1)/float(n-2) lam1=b0 lam2=2.0*b1-b0 lam3=6.0*b2-6.0*b1+b0 return lam1,lam2,lam3 def LN3(n,x,alpha): xm,xs,xr=MSR(n,x) b=1.0+xr*xr/2.0 t=(b+sqrt(b*b-1.0))**(1/3)+(b-sqrt(b*b-1.0))**(1/3)-1.0 ym=log(xs/sqrt(t*(t-1.0))) ys=sqrt(log(t)) aa=xm-exp(ym)*exp(ys*ys/2.0) pp=NEP(n,alpha) ww=np.array(norm.ppf(pp,0,1)) x_est=aa+np.exp(ym+ys*ww) return aa,ym,ys,pp,x_est def LP3(n,x,alpha): y=np.log(np.array(x)) ym,ys,yr=MSR(n,y) bb=4.0/yr**2 aa=ys/sqrt(bb) if yr<0.0: aa=-ys/sqrt(bb) cc=ym-aa*bb pp=NEP(n,alpha) ww=np.array(gamma.ppf(pp,bb)) x_est=np.exp(cc+aa*ww) return aa,bb,cc,pp,x_est def GEV(n,x,alpha): lam1,lam2,lam3=LMOM(n,x) d=2.0*lam2/(lam3+3.0*lam2)-log(2)/log(3) kk=7.8590*d+2.9554*d*d aa=kk*lam2/(1.0-2.0**(-kk))/Gamma(1.0+kk) cc=lam1-aa/kk*(1.0-Gamma(1.0+kk)) pp=NEP(n,alpha) x_est=cc+aa/kk*(1.0-(-np.log(pp))**kk) return kk,aa,cc,pp,x_est def DRAWFIG(x_org,x_est,str_t): xmin=0; xmax=400 ymin=0; ymax=400 fig=plt.figure() ax = fig.add_subplot(111) ax.plot(x_org,x_est,color='white',marker='o',markersize=7,markeredgecolor="b",markeredgewidth=1) ax.plot([xmin,xmax],[ymin,ymax],'-',color='black') ax.text(0.05,0.95,str_t,fontsize=12,fontweight='bold',color='black',ha='left',va='top',transform=ax.transAxes) txt='r={0:7.4f}'.format(np.corrcoef(x_org,x_est)[0,1]) ax.text(0.05,0.90,txt,fontsize=12,color='black',ha='left',va='top',transform=ax.transAxes) ax.set_xlim(xmin,xmax) ax.set_ylim(ymin,ymax) ax.set_xlabel('Observed data (mm/day)') ax.set_ylabel('Estimated data (mm/day)') aspect=1.0*(ymax-ymin)/(xmax-xmin)*(ax.get_xlim()[1]-ax.get_xlim()[0])/(ax.get_ylim()[1]-ax.get_ylim()[0]) ax.set_aspect(aspect) plt.savefig('fig_qq.png',dpi=200) plt.show() plt.close() # Main routine alpha=0.4 param=sys.argv fnameR=param[1] fnameW=param[2] # Input data x_org=[] fin=open(fnameR,'r') text=fin.readline() while 1: text=fin.readline() if not text: break text=text.strip() text=text.split() x_org=x_org+[float(text[0])] fin.close() n=len(x_org) x_org.sort() ik=0 if ik==0: str_t='LN3'; aa,ym,ys,pp,x_est=LN3(n,x_org,alpha) if ik==1: str_t='LP3'; aa,bb,cc,pp,x_est=LP3(n,x_org,alpha) if ik==2: str_t='GEV'; kk,aa,cc,pp,x_est=GEV(n,x_org,alpha) # Output data fout=open(fnameW,'w') print('{0:>5s} {1:>10s} {2:>10s} {3:>10s}'.format('i','p','observed','estimated'),file=fout) for i in range(0,n): print('{0:5d} {1:10.3f} {2:10.1f} {3:10.1f}'.format(i+1,pp[i],x_org[i],x_est[i]),file=fout) fout.close() DRAWFIG(x_org,x_est,str_t)