No.08917　極値統計(Gumbel分布)についての質問　　【HN】　2009/01/14(Wed) 15:06

現在仕事である製品(平板状）の上にある傷の最大の長さを評価ようとしてます。

計測器の関係で1度に測定できる面積は決まっています。
また，製品の全面を計測できれば問題は簡単なのですが，設備，工数の制約上現実的ではありません。

そこで極値統計(Gumbel分布)を利用して傷の最大値を推定しようと考えております。

手法は以下のとおり

(1)，検査基準面積(測定視野面積)を決定する。
(2)，計測器を使い検査基準面積内に存在する傷の最大の長さのものを選びその長さを計測する。
(3)，上記(2)を場所が重複しないようにN回繰り返し，N個のデータを抜き出す。
(この際，このN個のデータはGumbel分布に従う。)
(4)，(3)で得たデータ(長さ)を小さい順に並び替える。(L(1), L(2), L(3).....L(N))
(5)，基準化変数Yを以下のようにとる。
Y=-ln[-ln{i/(n+1)}]
i：(4)で並び替えたデータ(長さ)の小さいほうからの順番
n：全データ数

(6)，直行座標系の(X, Y)に(L(i), Y)をプロットし回帰直線を導出する。
(このデータの分布がGumbel分布に従ってるなら直線に近似できる。）

(7)，(6)で算出された回帰直線からある製品上の傷の最大長さを推定する。
（例：観測基準面積の1000倍の面積を持つ製品の最大傷長さは(6)の回帰直線でY=-ln[-ln{(1000/1001)}]=6.91になるXの値となる。）

以上のような手法である面積の最大傷長さを推定しようと思いますが，疑問があります。

面 積が1000である製品を検査基準面積1で50個のデータを取ったときにその製品の最大傷長さは回帰直線で Y=-ln[-ln{(1000/1001)}]=6.91でのXの値になりますが，同じ製品に検査基準面積2で50個のデータをとった場合，最大傷長さ は回帰直線でY=-ln[-ln{(500/501)}]=6.21でのXの値になってしまいます。

もし，データ数が十分大きいならば観測面積が1と2の場合のデータ分布は同一になっていきますが推定量を算出するためのYの値が異なってしまいます。

そこで質問です。
I，なぜこのような差が出るのか？(恐らく，ある分布から抜き出した最大値の分布はGumbel分布になるというところが理解できてないかと思います。）
II，この場合検査基準面積1と2の場合どちらが確かな値が出るのか？また，その判定基準は？。

近くの図書館でGumbel分布(極値統計)の参考書を探してみましたがこれといったものが見つからず途方にくれています。どなたかアドバイスいただけませんか？
よろしくお願いします。

以下，参考にした資料
【参考にした文献】装置材料の寿命予測入門 -極値統計の腐食への適用- 腐食防食協会編 丸善株式会社
【参考URL】http://www.bio.hi-tech.ac.jp/koyama/pdf/kyoku/kt2.htm　極値統計法による腐食材の残存寿命推定　【東電工業株式会社委託研究報告書,1994.3】 小山　信次
No.08936　Re: 極値統計(Gumbel分布)についての質問　　【Yoshi】　2009/01/17(Sat) 01:26

ここでコメントするのは初めてですが，誰も回答しないようですので答えようと思います．
まず，質問Iにつ いてですが，極値データというのはある期間（領域）の最大値をとることを繰り返したときに，得られるものなので，検査面積が異なれば再現確率（によって計 算される基準化変数）が変わるのは当然です．HNさんが対象としている全体で1000の面積の最大傷長さを調べるには，検査面積が1ならそれを1000回 検査したときに最大の長さとなったものは，1/1000の超過確率(非超過確率=1-1/1000≒1000/1001)で発生する最大傷長さということ になります．一方，検査面積が2なら，500回検査したときに最大の長さとなったものなので，500回に1度の超過確率(非超過確 率=1-1/500≒500/501)の最大傷長さになります．検査面積2の方が広いため，より長い最大傷長さで構成される極値データが集まるので，求ま るGumbel分布のパラメータは当然異なります．従って，検査面積1と同じ非超過確率で計算したのでは間違った再現確率の最大傷長さを計算してしまうこ とになります．
そして，質問IIの件ですが，検査するサンプル数が同じなのであれば，外挿によって求める再現確率が小さい方が精度は高くなると考 えられるので，（労力があまり変わらないのであれば）検査面積2でデータを取った方がいいと思います．このことは，Rを使うと次のように確認できます．
install.packages("evir")	#極値統計用のパッケージです
library(evir)
rgum <- function(n,mu=0,beta=1) -log(-log(runif(n)))*beta+mu	#Gumbel分布の逆関数から乱数を作る関数
rep1 <- rep2 <- rep3 <- rep(NA,1000)	#結果の入れ物
for(i in 1:1000){
	x <- rgum(100,2,0.5)	#Gumbel分布から発生させた100個の乱数
	x1 <- x[1:50]	#前半の50個を検査面積1のサンプルだとする
	x2 <- apply(cbind(x[1:50],x[51:100]),1,max)	#前半の50個と対応する後半の50個を比べて大きいものを選ぶ（検査面積2のサンプルだとする）
	par1 <-gumbel(x1)$par.ests	#最尤法による当てはめ（検査面積=1）
	par2 <-gumbel(x2)$par.ests	#最尤法による当てはめ（検査面積=2）
	F1 <- 1-1/1000	#1/1000の非超過確率（累積密度）（検査面積1）
	F2 <- 1-1/500	#1/500の非超過確率（累積密度）（検査面積2）
	rep1[i] <- -log(-log(F1))*par1[1]+par1[2]	#全体の最大傷長さの推定量（検査面積1）
	rep2[i] <- -log(-log(F2))*par2[1]+par2[2]	#全体の最大傷長さの推定量（検査面積2）
	rep3[i] <- -log(-log(F1))*par2[1]+par2[2]	#検査面積2なのに検査面積1の再現期間確率を用いた場合
}
mean(rep1)#検査面積1の乱数生成->推定を1000回繰り返した平均値
 [1] 5.3974
mean(rep2)#検査面積2の乱数生成->推定を1000回繰り返した平均値
 [1] 5.420649
mean(rep3)#検査面積2で違う再現期間超過確率を与えたときの平均値
 [1] 5.763179
-log(-log(1-1/1000))*0.5+2	#乱数を発生させたモデルの真の最大傷長さ
 [1] 5.453628
こ の場合は，検査面積2の方が真の最大傷長さに近い値となっていますが，検査面積1の方が真の最大傷長さに近くなる場合もあります．やはり，最も重要なのは サンプル数だと思います．そして，さらに重要なのはそのサンプルが本当にGumbel分布に合っているのかということだと思います．

私も 極値統計を勉強しているのですが，決まった参考書が無く苦労しております．しかし，結局は統計学の一部なので，「極値統計」にこだわらずに色々試してみる しかないと思います．また，確率紙による当てはめというのは若干古い手法のように思われます．Rを使えとは言いませんが，使えると上記のevirの gumbel関数のように，最尤法による当てはめなどが簡単に行えて便利ですので，一度勉強してみるといいと思います．水文統計でよければ以下の文献が参 考になるのではないかと思います．"年"を"面積"と置き換えて読んで下さい．
http://thesis.ceri.go.jp/center/doc/geppou/ceri/0005005050.pdf

● 「統計学関連なんでもあり」の過去ログ--- 042 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る