草薙の研究ログ

英語教育関係。でも最近は統計(特にR)ネタが中心。

一般化パレート分布をデータに当てはめる

一般化パレート分布は所得の分布などに使われるそうだ。
外国語教育研究でもこういった分布になる変数を私はひとつだけ知っている(いわないwww)。

Rにいろいろあると思うけど,ここではactuarパッケージとfitdistrplusパッケージを使う。
actuarパッケージに関数があるから,それをfitdistrplusパッケージのfitdist関数当てはめるというわけ。

#パッケージの準備
library(actuar)
library(fitdistrplus)

#乱数作っちゃう
#第一形状母数が3,第二形状母数が3,尺度母数が1000
set.seed(1)
dat<-rgenpareto(1000,3,3,scale=1000)

#経験分布を可視化
par(mfrow=c(1,2))
plot(ecdf(dat),main="")
hist(dat,col="lightblue",main="")


f:id:kusanagik:20170131170612p:plain

#最尤推定(初期値は適当)
fit<-fitdist(dat,"genpareto",start=list(shape1=1,shape2=1,scale=1000))
summary(fit)

Fitting of the distribution ' genpareto ' by maximum likelihood 
Parameters : 
          estimate  Std. Error
shape1    2.973575   0.3122382
shape2    2.811662   0.2800328
scale  1070.970526 236.9332116
Loglikelihood:  -8230.236   AIC:  16466.47   BIC:  16481.2 
Correlation matrix:
           shape1     shape2      scale
shape1  1.0000000 -0.6723317  0.9114599
shape2 -0.6723317  1.0000000 -0.9020353
scale   0.9114599 -0.9020353  1.0000000

plot(fit)


f:id:kusanagik:20170131170939p:plain


うむ。