evd 極値分布パッケージ

  • 繰り返し観察の極端な値(最大1時間降水量とか。多重検定の最大統計量とかも・・・)の確率密度分布は極値分布と呼ばれる分布になることが知られている
  • その形は、一般に、3つのパラメタ(location:\mu,scale:\sigma,shape:\zeta)で規定される、以下の式で表される
    • G(x)=exp^{-(1+\zeta(\frac{x-\mu}{\sigma}))^{-\frac{1}{\zeta}}}
  • この一般式の特別な場合として3つの名前のついた極値分布がある
    • Gumbel分布
    • Fure'chet分布
    • Weibull分布
  • Rではevdというパッケージがこの分布の関数を提供するとともに
    • dgev(x, loc=0, scale=1, shape=0, log = FALSE)
    • pgev(q, loc=0, scale=1, shape=0, lower.tail = TRUE)
    • qgev(p, loc=0, scale=1, shape=0, lower.tail = TRUE)
    • rgev(n, loc=0, scale=1, shape=0)
  • 観測データから、その極値分布パラメタを最尤法によって推定する関数fgevを提供している
  • 構造化による統計量分散インフレ(g:gamma)のあるようなカイ自乗統計量について、N回の検定を自由度dfにて行うことをiter回繰り返したときの、最大統計量の分布について、極値分布パラメタの推定を行い、その当てはまりの様子をプロットしてみるのが以下のソースで、掲載図(黒はシミュレーション。赤は推定極値分布からの乱数)
iter<-100
N<-10000
g<-2
df<-2
chis<-matrix(rchisq(iter*N,df)*g,nrow=iter)
plot(ppoints(N,a=0),sort(pchisq(chis[1,],df,lower.tail=TRUE)),type="l")
plot(log(ppoints(N,a=0),10),sort(log(pchisq(chis[1,],df,lower.tail=TRUE),10)),type="l")

maxchis<-apply(chis,1,max)
gevest<-fgev(maxchis)
gevest
rfromgev<-rgev(iter,loc=gevest$estimate[1],scale=gevest$estimate[2],shape=gevest$estimate[3])

plot(ppoints(iter,a=0),sort(maxchis),ylim=c(0,max(maxchis)),type="l")
par(new=T)
plot(ppoints(iter,a=0),sort(rfromgev),ylim=c(0,max(maxchis)),type="l",col="red")