2015-02-16

対数正規分布

対数正規分布 R

対数正規分布は、supportが正の値をとる連続分布で、生命現象を含む自然現象に広くみられる分布(こちらのFacultyOf1000 Prime)にもあるように
Rの対数正規分布に関する関数はdlnorm(),plnorm(),qlnorm(),rlnorm()の一式なのだが、ちょっと使いづらいかもしれない
どうして使いにくいか、というと、関数の引数が、「正規分布の平均がMで標準偏差がSであるときに、それを対数化した分布をdlnorm(x,log(M),log(S))として計算する」というような作りになっているから
これだと、「対数正規分布の平均がmeaで最頻値がmodであるような分布について調べたい」、という作業がしにくい
まず、「正規分布の平均がMで標準偏差がSであるときに、それを対数化した分布をdlnorm(x,log(M),log(S))として計算する」について確認しておく
- plnorm()はpnorm()と相性がよい

n <- 10^3
t <- runif(n,min=-10,max=10)
t.l <- exp(t)
M <- runif(length(t))
S <- runif(length(t))
dn <- pnorm(t,M,S)
dl <- plnorm(t.l,M,S)
plot(dn,dl)
abline(0,1,col=2)

- 正規分布は左右対称で対数正規分布は非対称であるから、pxxx()が対応づくということは、dxxx()は対応づかないということ
次に、「対数積分布の平均がmeaで最頻値がmodであるような分布」を扱うための関数を作っておく
- こちらにあるように、対数正規分布の形を決めるパラメタ(M,S)を用いて、対数正規分布の期待値mea、最頻値mod、中央値medは
- $mea = e^{M+\frac{S^2}{2}}$
- $mod = e^{M-S^2}$
- $med = e^M$
- となる。最後の中央値が $e^M$ となっていることが、「 $\log{mea}=M$ であるから、平均がMの正規分布を対数化したもの」ということを指していて、中央値・クオンタイルが正規分布と相性がよいことを表している
- medをmea,modで表せば $med = (mea^2\times mod)^{1/3}$
- M,Sをmea,modで表せば
  - $S=\sqrt{2/3 \log{mea/mod}}$
  - $M=\log{mod} + S^2$
- したがって

my.dlnorm <- function(x,mea=1,mod=0.5){
	S <- sqrt(2/3*log(mea/mod))
	M <- log(mod) + S^2
	dlnorm(x,M,S)
}
x <- seq(from=0,to=200,length=101)

mea <- 150
mod <- 100
d <- my.dlnorm(x,mea=mea,mod=mod)
plot(x,d,type="l")
abline(v=mea,col=2)
abline(v=mod,col=3)

my.plnorm <- function(q,mea=1,mod=0.5,lower.tail=TRUE){
	S <- sqrt(2/3*log(mea/mod))
	M <- log(mod) + S^2
	plnorm(q,M,S,lower.tail=lower.tail)
}

q <- x

p <- my.plnorm(q,mea=mea,mod=mod)
plot(q,p,type="l")
abline(v=mea,col=2)
abline(v=mod,col=3)

my.qlnorm <- function(p,mea=1,mod=0.5,lower.tail=TRUE){
	S <- sqrt(2/3*log(mea/mod))
	M <- log(mod) + S^2
	qlnorm(p,M,S,lower.tail=lower.tail)
}
p <- seq(from=0,to=1,length=101)

q <- my.qlnorm(p,mea=mea,mod=mod)
plot(p,q,type="l")
abline(h=mea,col=2)
abline(h=mod,col=3)

my.rlnorm <- function(n,mea=1,mod=0.5){
	S <- sqrt(2/3*log(mea/mod))
	M <- log(mod) + S^2
	rlnorm(n,M,S)
}

n <- 10^5
r <- my.rlnorm(n,mea=mea,mod=mod)
hist(r)
abline(v=mea,col=2)
abline(v=mod,col=3)

shinyでインターラクティブに

---
title: "log-normal distribution"
author: "ryamada"
date: "Monday, February 16, 2015"
output: html_document
runtime: shiny
---

```{r, echo=FALSE}
inputPanel(
  sliderInput("mod2", label = "mode:",
              min = 0, max = 500, value = 100, step = 1)
,
sliderInput("mean2", label = "mean:",
              min = 0, max = 1000, value = 200, step = 1)
,
sliderInput("xmax2", label = "xmax:",
              min = 0, max = 1000, value = 200, step = 1)
)

renderPlot({
xmax <- input$xmax2
mea <- input$mean2
mod <- input$mod2
med <- (mea^2*mod)^(1/3)

S <- sqrt(2/3*log(mea/mod))
M <- log(mod) + S^2

t <- seq(from=0,to=xmax,length=101)
d <- dlnorm(t,M,S)
plot(t,d,type="l")
abline(v=mea,col=2)
abline(v=mod,col=3)
abline(v=med,col=4)

})
```