草薙の研究ログ

英語の先生をやってます。

左に歪んだデータへのフィッティング

確率分布のフィッティングでは,右に歪んだデータを対象にすることが多いのだけど,たまに左に歪んだ分布が手に入ることもある(てか私は最近はじめて手にいれた)。いろいろあるんだろうけど,一般化極値分布でいいかな。えっと,こんなデータだとしよう。

f:id:kusanagik:20170711210454p:plain

#準備とデータの作成
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)

f:id:kusanagik:20170711210919p:plain

#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)

f:id:kusanagik:20170711211154p:plain

#無駄に無情報事前分布で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)

f:id:kusanagik:20170711211651p:plain