対数正規分布から、平均値(中央値)と標準偏差を指定して乱数を生成させる

R で対数正規分布から乱数を生成するには、rlnrom() 関数を使います。

data <- rlnorm(10000, meanlog = 0, sdlog = 1)
hist(data)

このとき、指定する対数正規分布のパラメータ meanlog, sdlog は、分布自身の平均値、標準偏差に対応していません。

cat("mean = ", mean(data), ", sd = ", sd(data), "\n", sep = "")

## mean = 1.651, sd = 2.093

シミュレーションなどで、対数正規分布自身の平均値、標準偏差を指定して乱数を発生させたいときがあります。

本記事では、対数正規分布自身の平均値(または中央値)、標準偏差を指定して乱数を発生させる方法について説明します。

1. 平均値、標準偏差を指定した乱数の生成

生成したい対数正規分布の平均値 \(E\)、標準偏差 \(S\) が与えられたとき、それに対する対数正規分布のパラメータ \(\mu\), \(\sigma\) を求めます。

それぞれの理論値は

\[ \begin{eqnarray} E &=& e^{\mu+\frac{\sigma^2}{2}} \\ S^2 &=& e^{2\mu+\sigma^2}(e^{\sigma^2}-1) \end{eqnarray} \]

となります。
まずは \(S\) の式を変形していきます。

\[ \begin{eqnarray} S^2 &=& e^{2\mu+\sigma^2}(e^{\sigma^2}-1) \\ &=& E^2(e^{\sigma^2}-1) \\ \frac{S^2}{E^2} &=& e^{\sigma^2} - 1 \\ e^{\sigma^2} &=& \frac{S^2}{E^2} + 1 \\ \sigma^2 &=& \log(\frac{S^2}{E^2} + 1) \\ \sigma &=& \sqrt{\log(\frac{S^2}{E^2} + 1)} \end{eqnarray} \]

\(E\) の式に対しては、

\[ \begin{eqnarray} E &=& e^{\mu+\frac{\sigma^2}{2}} \\ \log E &=& \mu + \frac{\sigma^2}{2} \\ \mu &=& \log E + \frac{\sigma^2}{2} \end{eqnarray} \]

最終的に

\[ \begin{eqnarray} \mu &=& \log E + \frac{\sigma^2}{2} \\ \sigma &=& \sqrt{\log\big(\frac{S^2}{E^2} + 1\big)} \end{eqnarray} \]

が得られました。

R で実行するには下記のようにします。

rlnorm2 <- function(n, mean=1, sd=1) {
  sdlog <- sqrt(log((sd/mean)^2 + 1))
  meanlog <- log(mean) - (sdlog^2) / 2
  rlnorm(n, meanlog = meanlog, sdlog = sdlog)
}
data <- rlnorm2(10000, mean=5, sd=2)
hist(data)

cat("mean = ", mean(data), ", sd = ", sd(data), "\n", sep = "")

## mean = 4.994, sd = 1.955

2. 中央値、標準偏差を指定した乱数の生成

次に、中央値を指定した場合について考えます。

生成したい対数正規分布の中央値 \(M\)、標準偏差 \(S\) が与えられたとき、それに対する対数正規分布のパラメータ \(\mu\), \(\sigma\) を求めます。

それぞれの理論値は

\[ \begin{eqnarray} M &=& e^{\mu} \\ S^2 &=& e^{2\mu+\sigma^2}(e^{\sigma^2}-1) \end{eqnarray} \]

となります。
まずは \(S\) の式を変形していきます。

\[ \begin{eqnarray} && S^2 = e^{2\mu+\sigma^2}(e^{\sigma^2}-1) \\ && S^2 = M^2 e^{\sigma^2}(e^{\sigma^2}-1) \\ && (e^{\sigma^2})^2 - e^{\sigma^2} - \frac{S^2}{M^2} = 0 \\ && e^{\sigma^2} = \frac{1 + \sqrt{4\frac{S^2}{M^2} + 1}}{2} \\ && \sigma^2 = \log \big( \frac{1 + \sqrt{4\frac{S^2}{M^2} + 1}}{2} \big) \\ && \sigma = \sqrt{ \log \big( \frac{1 + \sqrt{4\frac{S^2}{M^2} + 1}}{2} \big) } \end{eqnarray} \]

\(M\) については

\[ \begin{eqnarray} M &=& e^{\mu} \\ \log M &=& \mu \\ \mu &=& \log M \end{eqnarray} \]

最終的に

\[ \begin{eqnarray} \mu &=& \log M \\ \sigma &=& \sqrt{ \log \big( \frac{1 + \sqrt{4\frac{S^2}{M^2} + 1}}{2} \big) } \end{eqnarray} \]

が得られました。

R で実行するには下記のようにします。

rlnorm3 <- function(n, median=1, sd=1) {
  meanlog <- log(median)
  sdlog <- sqrt(log((1 + sqrt((2*sd/median)^2 + 1))/2))
  rlnorm(n, meanlog = meanlog, sdlog = sdlog)
}
data <- rlnorm3(10000, median=5, sd=2)
hist(data)

cat("median = ", median(data), ", sd = ", sd(data), "\n", sep = "")

## median = 4.995, sd = 2.047

3. まとめ

対数正規分布に対して、平均値(または中央値)、標準偏差を指定して乱数を生成する方法について説明しました。

4. 参考

対数正規分布 - Wikipedia