読者です 読者をやめる 読者になる 読者になる

草薙の研究ログ

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

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

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)


ちなみに,retimesパッケージにはtimefitっていう専用の関数がある。ま,ほぼおなじ推定値が帰ってくる。

timefit(dat)


なんかもうoptimやりゃいいって気がしてくるな。