左に歪んだデータへのフィッティング
確率分布のフィッティングでは,右に歪んだデータを対象にすることが多いのだけど,たまに左に歪んだ分布が手に入ることもある(てか私は最近はじめて手にいれた)。いろいろあるんだろうけど,一般化極値分布でいいかな。えっと,こんなデータだとしよう。
#準備とデータの作成 library(evd) library(fitdistrplus) library(ismev) library(MCMCpack) dat<-rgev(1000,20,3,-.5) hist(dat,main="",xlab="Value",breaks=20,col="pink") #まずはfitdist関数でフィット #初期値は適当に標本平均と標準偏差いれちゃう fit<-fitdist(dat,"gev",start=c(mean(dat),sd(dat),0)) summary(fit) Fitting of the distribution ' gev ' by maximum likelihood Parameters : estimate Std. Error 1 19.9521276 0.10240292 2 3.0110352 0.07695938 3 -0.4899288 0.01781197 Loglikelihood: -2398.873 AIC: 4803.746 BIC: 4818.469 Correlation matrix: [,1] [,2] [,3] [1,] 1.0000000 -0.2868652 -0.3588115 [2,] -0.2868652 1.0000000 -0.6712985 [3,] -0.3588115 -0.6712985 1.0000000 plot(fit)
#ismevパッケージでフィット fit2<-gev.fit(dat) fit2 $conv [1] 0 $nllh [1] 2398.873 $mle [1] 19.9520718 3.0115839 -0.4900573 $se [1] 0.10241636 0.07699955 0.01780562 gev.diag(fit.gev)
#無駄に無情報事前分布でMCMCしてみる llf<-function(beta,x){ sum(log(dgev(x,beta[1],beta[2],beta[3]))) } m<-MCMCmetrop1R(llf,theta.init=c(mean(dat),sd(dat),0),x=dat) summary(m) Iterations = 501:20500 Thinning interval = 1 Number of chains = 1 Sample size per chain = 20000 1. Empirical mean and standard deviation for each variable, plus standard error of the mean: Mean SD Naive SE Time-series SE [1,] 19.9413 0.10164 0.0007187 0.0024210 [2,] 3.0184 0.07621 0.0005389 0.0018702 [3,] -0.4871 0.01757 0.0001242 0.0004275 2. Quantiles for each variable: 2.5% 25% 50% 75% 97.5% var1 19.7453 19.871 19.9407 20.0099 20.1407 var2 2.8700 2.966 3.0172 3.0689 3.1709 var3 -0.5214 -0.499 -0.4873 -0.4754 -0.4516 plot(m)