草薙の研究ログ

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

ジョンソンのSU分布

ジョンソンの分布っていうのは,こんな感じ。PDFとかCDFとか。

Johnson's SU-distribution - Wikipedia, the free encyclopedia


この分布は4パラミターある。

  1. γ
  2. δ
  3. ξ
  4. λ

このように贅沢な関数なので,やはり表現力は高い。
Rでは,SuppDistsパッケージなどでこれを扱える。

https://cran.r-project.org/web/packages/SuppDists/SuppDists.pdf


たとえば,微妙に歪んだこんな分布を考える(ex-Gaussian分布から乱数生成)。

#psychパッケージとretimesパッケージをつかう
library(psych)
library(retimes)

#乱数を生成
d<-rexgauss(1000,3000,200,1000)

#モーメントを確認
describe(d)

#ヒストグラムで可視化
hist(d,col="lightgray")

f:id:kusanagik:20160826182715p:plain

#ジョンソン分布へフィットさせてパラミタを得る
param<-JohnsonFit(d)

#パラミタを確認
param

#このパラミタにもとづいて乱数を生成
s<-rJohnson(1000,param)

#この乱数をヒストグラムで可視化
hist(s,col="lightgray")

f:id:kusanagik:20160826183141p:plain

ここからex-Gaussianのパラミタを得てみる。

timefit(s)

私の場合だと,μ = 2998,σ = 219,τ = 976で,正解がそれぞれ3000, 200, 1000だから,まあいい感じだった。

確率密度関数を求めるには,dJohnsonっていう関数を使う。

x<-seq(0,10000,.01)
plot(x,dJohnson(x,param),type="l",ylab="Probability",ylim=c(0,.001))
lines(x,dexgauss(x,3000,200,1000),type="l",ylab="Probability",col=2)

f:id:kusanagik:20160826183911p:plain

赤はex-Gaussian分布の関数。

いい感じ。