#********************************** # Probability paper #********************************** BEGIN{ ind=ARGV[1] FILENAME=ARGV[2] delete ARGV[1] } { if(2<=NR)data[NR-1]=$1 } END{ nd=NR-1 #-------------------------- # Sort (small->large) #-------------------------- for(i=1;i<=nd;i++){ for(j=i+1;j<=nd;j++){ if(data[i]>data[j]){ work=data[i];data[i]=data[j];data[j]=work; } } } #-------------------------- # Non-exceeding probability # Plotting and positioned formula # Weibull : alpha=0 # Cunnane : alpha=0.4 # Hazen : alpha=0.5 #-------------------------- alpha=0.4 for(i=1;i<=nd;i++){ p=(i-alpha)/(nd+1-2*alpha) if(ind==1){datax[i]=data[i] ;datay[i]=TODA(p, ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);} if(ind==2){datax[i]=log(data[i])/log(10);datay[i]=TODA(p, ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);} if(ind==3){datax[i]=data[i] ;datay[i]=GUMBEL(p);} if(ind==4){datax[i]=log(data[i])/log(10);datay[i]=WEIBLL(p);} if(ind==5){datax[i]=data[i] ;datay[i]=EXPO(p);} } #-------------------------- # Regression (y=aa*x+bb) #-------------------------- x1=0.0;y1=0.0;x2=0.0;xy=0.0; for(i=1;i<=nd;i++){ x1=x1+datax[i]; y1=y1+datay[i]; x2=x2+datax[i] * datax[i]; xy=xy+datax[i] * datay[i]; } xm=x1/nd;ym=y1/nd; aa=(nd*xy-x1*y1)/(nd*x2-x1*x1); bb=(x2*y1-x1*xy)/(nd*x2-x1*x1); c1=0.0;c2=0.0;c3=0.0; for(i=1;i<=nd;i++){ c1=c1+(datax[i]-xm)*(datay[i]-ym); c2=c2+(datax[i]-xm)*(datax[i]-xm); c3=c3+(datay[i]-ym)*(datay[i]-ym); } rr=c1/sqrt(c2*c3); #-------------------------- # Print out #-------------------------- print "#y=a*x+b" printf "#a=%4.3e\n",aa printf "#b=%4.3e\n",bb printf "#r=%4.3f\n",rr printf "#nd=%d\n",nd print "#*** Regression line ***" print "#",data[1],aa*datax[1]+bb print "#",data[nd],aa*datax[nd]+bb print "#*** y-axis data (text) ***" if(ind<=4){ prb[1]=0.999; prb[2]=0.99; prb[3]=0.90; prb[4]=0.80; prb[5]=0.70; prb[6]=0.60; prb[7]=0.50; prb[8]=0.40; prb[9]=0.30; prb[10]=0.20; prb[11]=0.10; prb[12]=0.01; prb[13]=0.001; } if(ind==5){ prb[1]=0.999; prb[2]=0.99; prb[3]=0.90; prb[4]=0.80; prb[5]=0.70; prb[6]=0.60; prb[7]=0.50; prb[8]=0.40; prb[9]=0.30; prb[10]=0.20; prb[11]=0.001; prb[12]=0.001; prb[13]=0.001; } nj=13 for(i=1;i<=nj;i++){ p=prb[i] if(ind==1)pY[i]=TODA(p, ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10) if(ind==2)pY[i]=TODA(p, ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10) if(ind==3)pY[i]=GUMBEL(p) if(ind==4)pY[i]=WEIBLL(p) if(ind==5)pY[i]=EXPO(p) } for(i=1;i<=nj;i++){ print "#",-0.02,pY[i],14,0,0,"MR",prb[i] } print "#*** y-axis data (grid) ***" for(i=1;i<=nj;i++){ print "#>" printf "#0 %g\n",pY[i] printf "#1 %g\n",pY[i] } print "#*** Measured data ***" for(i=1;i<=nd;i++){ printf "%g %g\n",data[i],datay[i] } } #################################################################### function TODA(rx, ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10) { # -------------------- # Toda's approxmation # -------------------- if(rx==0.5){ uay=0.0; }else{ if(rx<0.5){ ry=rx; }else{ ry=1-rx; } ay=-log(4.0*ry*(1.0-ry)); b0 = 1.570796288; b1 = 0.03706987906; b2 =-0.0008364353589; b3 =-0.0002250947176; b4 = 0.000006841218299; b5 = 0.000005824238515; b6 =-0.00000104527497; b7 = 0.00000008360937017; b8 =-0.000000003231081277; b9 = 0.00000000003657763036; b10= 0.0000000000006936233982; uay=sqrt(ay*(b0+b1*ay+b2*ay^2+b3*ay^3+b4*ay^4+b5*ay^5+b6*ay^6+b7*ay^7+b8*ay^8+b9*ay^9+b10*ay^10)); if(rx<0.5)uay=-uay; } return uay; } #################################################################### function GUMBEL(p) { return -log(-log(p)) } #################################################################### function WEIBLL(p) { return log(-log(1-p)) } #################################################################### function EXPO(p) { return -log(1-p) } ####################################################################