Rで逆正規分布の確率密度関数，累積分布関数，乱数発生

逆正規分布（逆ガウス分布，Inverse Gaussian Distribution）は，反応時間の分析で使ったりすることがあるらしい。歪んでいるので。

2母数の連続型分布で，形状パラミタ（λ）と平均（μ）をもつ。
相変わらず英語のwiki先生は統計と数学に詳しい。

Inverse Gaussian distribution - Wikipedia, the free encyclopedia

Rで関数関係がちょっと見て見つからなかったので書きました。

#確率密度関数
dig<-function(x,mu,rambda){
	(rambda/(2*pi*(x^3)))^(1/2)*exp((-1*rambda*((x-mu)^2))/(2*(mu^2)*x))
	}

#累積分布関数
pig<-function(x,mu,rambda){	
	pnorm(sqrt(rambda/x)*(x/mu-1))exp(2*rambda/mu)*pnorm(-1*sqrt(rambda/x)*(x/mu+1))
	}	

#乱数発生
rig<-function(n,mu,rambda){
	y<-rnorm(n)^2
	x<-mu+((mu^2)*y)/(2*(rambda))-(mu/(2*rambda))*sqrt((4*mu*rambda*y)+(mu^2)*(y^2))
	z<-runif(n,0,1)
	ifelse(z<=mu/(mu+x),x,(mu^2)/x)
}

#乱数発生のアルゴリズムはhttps://en.wikipedia.org/wiki/Inverse_Gaussian_distributionを参照
#誤りがあったらこっそり教えてください

…ってやっていたら，Rのstatmodっていうパッケージがあるのを発見した。

Gordon Smyth, Yifang Hu, Peter Dunn, Belinda Phipson and Yunshun Chen (2015). statmod: Statistical Modeling. R package version 1.4.22. http://CRAN.R-project.org/package=statmod

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

なので，興味がある方はこっちを使ってください。チェックがてら使ってみよう。

#呼び出し
library(statmod)

#確率密度関数
dig(1,1,1)	#私の関数
dinvgauss(1,mean=1,shape=1)	#statmodの関数

#累積分布関数
pig(1,1,1)
pinvgauss(1,mean=1,shape=1)


#乱数生成
r.mine<-rig(10000,1,1)
r.statmod<-rinvgauss(10000,mean=1,shape=1)

#10分位点で大雑把に比較
quantile(r.mine,seq(0,1,.1))
quantile(r.statmod,seq(0,1,.1))

#ヒストグラムで見比べ
par(mfrow=c(1,2))	#描画画面の分割
hist(r.mine)
hist(r.statmod)