ReturnPeriod(再現期間)

　　Return Period (再現期間) 河川屋 2006/07/14 23:42:14 
　　└つづき。 河川屋 2006/07/14 23:44:35

　　Return Period (再現期間) 河川屋 2006/07/14 23:42:14  ツリーへ

Return Period (再現期間) 返事を書く

河川屋 2006/07/14 23:42:14

たまには、何となく実務計算らしいものを。
過去のデータを用い、今年の大雨は何年に1度クラスの雨だったかを計算するプログラムです。
よくみると、偏差値を計算し、偏差値から順位を推測
するのと全く同じ方法だったりします。

DECLARE EXTERNAL FUNCTION NORM.NORMSDIST
DECLARE EXTERNAL FUNCTION NORM.NORMSINV

DIM Y(100)
OPEN #1: NAME "RAIN.TXT" !各年最大日雨量
! データの最後にマイナスの値を入れておく。
LET N=0
INPUT #1: R
DO WHILE R>0
LET N=N+1
LET Y(N)=LOG10(R) !年最大値(=極値)は対数をとると正規分布になる。
INPUT #1:R
LOOP
CLOSE #1
PRINT N
CALL STDEVA(Y,N,MEAN,SIGMA)
PRINT "今年の大雨の雨量＝";
INPUT PROMPT "今年の大雨の雨量＝":R0
LET R0=LOG10(R0)
LET HENSA=(R0-MEAN)/SIGMA
LET PROB=1/(1-NORMSDIST(HENSA))
LET KETA=1-LOG10(PROB)
LET PROB=ROUND(PROB,KETA) !計算結果の有効桁は１桁強のため、それなりに丸める。
PRINT "今年の大雨は " ,PROB,"年に１回の大雨です"
STOP

SUB STDEVA(Y(),N,MEAN,SIGMA) !不偏標準偏差と平均を求める。
LET TX=0
LET TXX=0
FOR I=1 TO N
LET TX=TX+Y(I)
LET TXX=TXX+Y(I)*Y(I)
NEXT I
LET MEAN=TX/N
LET SIGMA=SQR((N*TXX-TX*TX)/N/(N-1))
END SUB
END

MODULE NORM !標準正規累積分布関数の近似計算
!MS-EXCELの同名関数と同じ動作ですが少し計算精度が悪い。
PUBLIC FUNCTION NORMSDIST
PUBLIC FUNCTION NORMSINV

EXTERNAL FUNCTION NORMSDIST(X)
LET SPI=SQR(2*PI)
LET W=1/(1+0.2316419*X)
LET D=EXP(X^2/(-2))/SPI
LET NORMSDIST=1.0-D*(0.3193815*W-0.3565638*W^2+1.781478*W^3-1.821256*W^4+1.330274*W^5)
END FUNCTION

EXTERNAL FUNCTION NORMSINV (P)
!Hastingsの式で第1近似(=ANS0)を求め、NEWTON法で精度を上げる。
LET PP=P
IF PP > 0.5 THEN LET PP=1-P
LET Z=SQR(-2.0*LOG(PP))
LET ANS0=Z-(0.27061*Z+2.30753)/((0.04481*Z+0.99229)*Z+1.0)
LET P0=1-NORMSDIST(ANS0)
LET SPI=SQR(2*PI)
LET D=EXP(ANS0^2/(-2))/SPI
LET ANS=ANS0+(P0-PP)/D
IF P >= 0.5 THEN LET NORMSINV=ANS ELSE LET NORMSINV=-ANS
END FUNCTION
END MODULE

　　└つづき。 河川屋 2006/07/14 23:44:35  ツリーへ

Re: Return Period (再現期間) 返事を書く

河川屋 2006/07/14 23:44:35

つづき。

標準偏差と平均を求める方法ですが、前投稿は説明用で、こちらが実務用です。
(前のでも間違いとはいえないけど....)
統計の本を見る場合、「確率紙」に言及している本なら、こういうやりかたも
載っています。(注:おなじみの方法と、標準偏差の推定値は一致しない。)

SUB STDEVA(Y(),N,MEAN,SIGMA) !不偏標準偏差と平均を求める(順序確率による。)

LET M=N !全年使用
! LET M=5 !上位5年使用 (裾の部分のみを合わせたいとき)
! LET M=INT(N/10+0.5) !上位1/10個を使用 (　〃　)

! 降順に並べ換える
FOR I=1 TO N-1
LET K=I
FOR J=I+1 TO N
IF Y(K)<Y(J) THEN LET K=J
NEXT J
IF K<>I THEN SWAP Y(I),Y(K)
NEXT I

! 最小２乗直線回帰により傾き(=SIGMA)、切片(=MEAN)を求める
LET TX=0
LET TXX=0
LET TXY=0
FOR I=1 TO M
LET PL=I/(N+1) !プロットポジション＝トーマス法
!LET PL=(2*I-1)/(2*N) !プロットポジション＝ヘイズン法
LET X=-NORMSINV(PL)
LET TX=TX+X
LET TY=TY+Y(I)
LET TXX=TXX+X*X
LET TXY=TXY+X*Y(I)
NEXT I
LET MEAN=(TXX*TY-TX*TXY)/(M*TXX-TX*TX)
LET SIGMA=(M*TXY-TX*TY)/(M*TXX-TX*TX)
END SUB

インデックスへ EXIT
新規発言を反映させるにはブラウザの更新ボタンを押してください。

Return Period (再現期間)	返事を書く
河川屋 2006/07/14 23:42:14
たまには、何となく実務計算らしいものを。過去のデータを用い、今年の大雨は何年に1度クラスの雨だったかを計算するプログラムです。よくみると、偏差値を計算し、偏差値から順位を推測するのと全く同じ方法だったりします。 DECLARE EXTERNAL FUNCTION NORM.NORMSDIST DECLARE EXTERNAL FUNCTION NORM.NORMSINV DIM Y(100) OPEN #1: NAME "RAIN.TXT" !各年最大日雨量 ! データの最後にマイナスの値を入れておく。 LET N=0 INPUT #1: R DO WHILE R>0 LET N=N+1 LET Y(N)=LOG10(R) !年最大値(=極値)は対数をとると正規分布になる。 INPUT #1:R LOOP CLOSE #1 PRINT N CALL STDEVA(Y,N,MEAN,SIGMA) PRINT "今年の大雨の雨量＝"; INPUT PROMPT "今年の大雨の雨量＝":R0 LET R0=LOG10(R0) LET HENSA=(R0-MEAN)/SIGMA LET PROB=1/(1-NORMSDIST(HENSA)) LET KETA=1-LOG10(PROB) LET PROB=ROUND(PROB,KETA) !計算結果の有効桁は１桁強のため、それなりに丸める。 PRINT "今年の大雨は " ,PROB,"年に１回の大雨です" STOP SUB STDEVA(Y(),N,MEAN,SIGMA) !不偏標準偏差と平均を求める。 LET TX=0 LET TXX=0 FOR I=1 TO N LET TX=TX+Y(I) LET TXX=TXX+Y(I)Y(I) NEXT I LET MEAN=TX/N LET SIGMA=SQR((NTXX-TXTX)/N/(N-1)) END SUB END MODULE NORM !標準正規累積分布関数の近似計算 !MS-EXCELの同名関数と同じ動作ですが少し計算精度が悪い。 PUBLIC FUNCTION NORMSDIST PUBLIC FUNCTION NORMSINV EXTERNAL FUNCTION NORMSDIST(X) LET SPI=SQR(2PI) LET W=1/(1+0.2316419X) LET D=EXP(X^2/(-2))/SPI LET NORMSDIST=1.0-D(0.3193815W-0.3565638W^2+1.781478W^3-1.821256W^4+1.330274W^5) END FUNCTION EXTERNAL FUNCTION NORMSINV (P) !Hastingsの式で第1近似(=ANS0)を求め、NEWTON法で精度を上げる。 LET PP=P IF PP > 0.5 THEN LET PP=1-P LET Z=SQR(-2.0LOG(PP)) LET ANS0=Z-(0.27061Z+2.30753)/((0.04481Z+0.99229)Z+1.0) LET P0=1-NORMSDIST(ANS0) LET SPI=SQR(2PI) LET D=EXP(ANS0^2/(-2))/SPI LET ANS=ANS0+(P0-PP)/D IF P >= 0.5 THEN LET NORMSINV=ANS ELSE LET NORMSINV=-ANS END FUNCTION END MODULE

Re: Return Period (再現期間)	返事を書く
河川屋 2006/07/14 23:44:35
つづき。標準偏差と平均を求める方法ですが、前投稿は説明用で、こちらが実務用です。 (前のでも間違いとはいえないけど....) 統計の本を見る場合、「確率紙」に言及している本なら、こういうやりかたも載っています。(注:おなじみの方法と、標準偏差の推定値は一致しない。) SUB STDEVA(Y(),N,MEAN,SIGMA) !不偏標準偏差と平均を求める(順序確率による。) LET M=N !全年使用 ! LET M=5 !上位5年使用 (裾の部分のみを合わせたいとき) ! LET M=INT(N/10+0.5) !上位1/10個を使用 (　〃　) ! 降順に並べ換える FOR I=1 TO N-1 LET K=I FOR J=I+1 TO N IF Y(K)<Y(J) THEN LET K=J NEXT J IF K<>I THEN SWAP Y(I),Y(K) NEXT I ! 最小２乗直線回帰により傾き(=SIGMA)、切片(=MEAN)を求める LET TX=0 LET TXX=0 LET TXY=0 FOR I=1 TO M LET PL=I/(N+1) !プロットポジション＝トーマス法 !LET PL=(2I-1)/(2N) !プロットポジション＝ヘイズン法 LET X=-NORMSINV(PL) LET TX=TX+X LET TY=TY+Y(I) LET TXX=TXX+XX LET TXY=TXY+XY(I) NEXT I LET MEAN=(TXXTY-TXTXY)/(MTXX-TXTX) LET SIGMA=(MTXY-TXTY)/(MTXX-TXTX) END SUB