自分で用意した分布をデータへフィットさせる

fitdistrplusパッケージのfitdist関数は，dnorm, pnormのように，dとpの関数が定義されていればどんなものでも指定することができる。自作でもいいってこと。ま，でも殆どの場合，自分程度が思いつくような分布はすでに用意されているってのがR。

たとえば，レイリー分布の場合。レイリー分布はVGAMパッケージに含まれている。レイリー分布のパラミタはひとつ（形状母数）。

#パッケージの読み込み
library(VGAM)

#例の作成
dat<-rrayleigh(1000,300)

#fitdist関数によるフィット
#初期値を指定してやらなければならない。とりあえずここではデータの中央値を入れてやる。
fit<-fitdist(dat,"rayleigh",start=list(scale=median(dat)))
summary(fit)

Fitting of the distribution ' rayleigh ' by maximum likelihood 
Parameters : 
      estimate Std. Error
scale 296.7956   4.692907
Loglikelihood:  -6652.58   AIC:  13307.16   BIC:  13312.07 

plot(fit)

f:id:kusanagik:20170117133218p:plain

よし。

じゃ，次は指数正規合成分布。これはretimesパッケージにある。ただし，dexgaussとrexgaussしかないので，pexgaussを自分で作ってしまう。

#パッケージ読み込み
library(retimes)

#pexgaussを作成
pexgauss<-function(x,mu,sigma,tau){
tau<-1/tau
u<-tau*(x-mu)
v<-tau*sigma
phi1<-pnorm(u,0,v)
phi2<-pnorm(u,v^2,v)
e<-exp(1)
p<-phi1-e^(-u+(v^2)/2+log(phi2))
p
}

#数値例の作成
dat<-rexgauss(1000,1000,500,1000)

#fitdist関数によるフィット
#初期値を指定してやらなければならない。とりあえずここではデータの正解を入れてやる。
fit<-fitdist(dat,"exgauss",start=list(mu=1000,sigma=500,tau=1000))
summary(fit)