- 昨日、対数正規分布についてメモした
- その離散版。ポアソン対数正規分布
- こちらのqRT-PCRとか、こちらの細菌検出とかで使われている
- ポアソン分布のパラメタ
が対数正規分布に従っているとして計算しているらしい
- poilogなるパッケージがあって、dpoilog(),rpoilog()という関数がある
- このdpoilog()のソースをこちらなどで確認すると、対数正規分布の確率密度分布を出して、その重みづけ平均としてポアソン対数正規分布の密度を計算しているようだ(以下のように、一致する)
library(poilog)
my.dpslnorm2 <- function(x,mean=0,sd=1,k=1){
z <- seq(from=0,to=50,length=10^k+1)
d <- dlnorm(z,mean,sd)
Lambda <- z
D <- dpois(x,Lambda)*(z[2]-z[1])
sum(D*d)
}
xs <- 0:10
ret <- ret1 <- rep(0,length(xs))
for(i in 1:length(xs)){
x <- xs[i]
mea <- 3
sd <- 2
ret[i] <- dpoilog(x,mea,sd)
ret２[i] <- my.dpslnorm２(x,mean=mea,sd=sd)
}
ret
ret２
> ret
[1] 0.06086260 0.05421382 0.04571094 0.03902130 0.03388362 0.02986360
[7] 0.02664534 0.02401495 0.02182648 0.01997801 0.01839654
> ret1
[1] 0.06086266 0.05421384 0.04571095 0.03902131 0.03388363 0.02986361
[7] 0.02664535 0.02401495 0.02182648 0.01997801 0.01839654
-
- rpoilog()は0という観測をデフォルトでは「返さない」ので、ちょっと気持ち悪い
> rpoilog(10,0,1)
[1] 3 1 1 2 1 4
> rpoilog(10,0,1,keep0=TRUE)
[1] 0 0 0 3 3 0 4 2 1 0
- ついでに、ＲにはC実装の離散積分関数Rdqagsというのがあるらしいことが、poilogのコードからわかりました(こちら)
- さて。昨日と同じで、対数正規分布のパラメタ指定は、対数正規分布の期待値となっていなくて、「元となる」正規分布のそれなので面倒くさい
- 対数正規分布の期待値と「最頻値」とを基準にポアソン対数正規分布のd/r関数を作っておこう
- ポアソン分布の期待値はそのパラメタλに等しいから、要するに対数正規分布の期待値をいくつにするか、とポアソン対数正規分布の期待値をいくつにするかは同じ
- まず、rpoilog()関数から乱数の期待値を確認しておこう
mea <- 0.5
sd <- 1
n <- 10^6
rr <- rpoilog(n,mea,sd,keep0=TRUE)
mean(rr)
exp(mea+sd^2/2)
> mean(rr)
[1] 2.726748
> exp(mea+sd^2/2)
[1] 2.718282
-
- 対数正規分布の期待値と最頻値を引数とする関数を作る
my.dpoislnorm <- function(x,mea=1,mod=0.5){
S <- sqrt(2/3*log(mea/mod))
M <- log(mod) + S^2
dpoilog(x,M,S)
}
my.rpoislnorm <- function(n,mea=1,mod=0.5){
S <- sqrt(2/3*log(mea/mod))
M <- log(mod) + S^2
rpoilog(n,M,S,keep0=TRUE)
}
n <- 10^4
Mea <- 50
Mod <- 30
RR <- my.rpoislnorm(n,Mea,Mod)
mean(RR)
tab.cnt <- tabulate(RR+1)
plot(0:max(RR),tab.cnt,type="l")
abline(v=Mod,col=2)
