【Rで非ガウス統計３】相乗対数正規モデルのパラメタ推定：モーメント法

乱流、株式指標、心拍ゆらぎなどの時系列では、その確率分布が正規分布に比べて裾の厚い非ガウス性を示すことがよくあります。今回は、そのような非ガウス分布を記述するモデルとして、Castaing が発達乱流の記述のために提案した相乗対数正規モデルを導入し、非ガウス性の強さを特徴づけるパラメタの推定法を紹介します。この相乗対数正規モデルは、

局所的にはガウス分布に従い、
局所標準偏差が対数正規分布に従ってランダムにゆらぐ

という確率過程を考えています。この過程は、確率変数 $X$ が、正規分布に従う２つの確率変数 $W$ と $Y$ を使って、　

$\displaystyle{ X = W e^{Y} }$

と表すことができます。ここで

$W$ ：平均 0 のガウス過程
$Y$ ： $W$ と独立なガウス過程で、分散が $\lambda^ 2$ 、平均が $- \lambda^ 2$

です。 $X$ の非ガウス性の程度は、 $\lambda^ 2$ によって完全に特徴づけられます。 $Y$ の平均を $- \lambda^ 2$ とするのは、

$\displaystyle{ \mathrm{Var} (X) = \mathrm{Var} (W) }$

を成り立たせるためで、 $e^ {Y}$ が非ガウス性のみに影響を与えるようにする工夫です。
　このモデルは、下の図に描かれているように、

$\lambda^ 2 = 0$ のとき、 $X$ がガウス分布、
$\displaystyle{\lambda^ 2 \gt 0}$ で、 $\displaystyle{\lambda^ 2}$ が大きいほど、裾が厚く、中心が尖った非ガウス分布、

になります。

相乗対数正規モデルの確率密度関数（赤実線）。破線は正規分布（ガウス分布）。右側のグラフは縦軸を対数スケールにとったもの。

1. モーメント法を用いたパラメタ推定
2. 推定量の使い方
- ■ Rで実装
3. 数値実験
4. まとめ
5. 付録：確率密度関数を描くRスクリプト
- 乗法log-normalモデルの確率密度関数の描画

1. モーメント法を用いたパラメタ推定

　ここでは、平均の違いや分散の違いではなく、非ガウス性という形状の違いのみを調べたいので、データの解析では、平均 0、分散 1 に標準化した確率密度関数（probability density function: PDF）を考えたほうが便利です。この場合、平均 0、分散 1 の相乗対数正規モデルの分布は、

$\displaystyle{ f_{\lambda}(x)=\int_{0}^{\infty}\frac{1}{\sqrt{2\pi\,\lambda}}\exp\!\left(-\frac{(\ln\sigma+\lambda^{2})^{2}}{2\lambda^{2}}\right)\frac{1}{\sqrt{2\pi}\,\sigma}\exp\!\left(-\frac{x^{2}}{2\sigma^{2}}\right)\,d(\ln\sigma) }$

で与えられます。
　このとき、この分布に従う確率変数 $X$ の絶対値の $q$ 次モーメントは、次のようになります：

$\displaystyle{ \langle |X|^q \rangle = 2^{q/2}\, \frac{\Gamma\!\left(\frac{q+1}{2}\right)}{\sqrt{\pi}}\, \exp\!\left(\frac{q(q-2)}{2}\lambda^2\right) \qquad (q>-1) }$

ここで、 $\langle \, \cdot \, \rangle$ は期待値を表します。
　この式を $\lambda^ 2$ について解けば、観測時系列の絶対値モーメント $\langle|x|^ q\rangle$ から $\lambda^ 2$ を推定するための式

$\displaystyle{ \widehat{\lambda}^2_q = \frac{2}{q(q-2)} \left[ \ln\!\left( \frac{\sqrt{\pi}\,\langle |x|^q\rangle}{2^{q/2}} \right)- \ln \Gamma \!\left(\frac{q+1}{2}\right) \right] }$

が導けます。ただし、この式では、 $q\neq 0, 2$ ですので、 $q = 0$ 、 $q = 2$ のときには、それぞれ、

$\displaystyle{ \begin{aligned} & \hat{\lambda}_0^2=-\langle\ln | x| \rangle-\frac{\gamma+\ln 2}{2} \quad(q=0) \\ & \hat{\lambda}_2^2=-\left\langle x^2 \ln \right| x| \rangle-\frac{\gamma+\ln 2}{2}-1 \quad(q=2) \end{aligned} }$

を使ってください。

2. 推定量の使い方

　先ほど導いた非ガウスパラメタの推定量

$\displaystyle{ \begin{aligned} \hat{\lambda}^2_q &=\frac{2}{q(q-2)} \left[ \ln\!\left( \frac{\sqrt{\pi}\,\langle |x|^q\rangle}{2^{q/2}} \right)- \ln \Gamma \!\left(\frac{q+1}{2}\right) \right] & (q \neq 0, 2) \\ \hat{\lambda}_0^2 &=-\langle\ln | x| \rangle-\frac{\gamma+\ln 2}{2} &(q=0) \\ \hat{\lambda}_2^2 &=-\left\langle x^2 \ln \right| x| \rangle-\frac{\gamma+\ln 2}{2}-1 &(q=2) \end{aligned} }$

いずれか一つの方法を用いれば、 $\lambda^ 2$ を推定することは可能ですが、このモデルを実際のデータに適用する際には、観測データが相乗対数正規モデルによって十分に近似できているかどうかを、必ず確認する必要があります。その確認方法としては、主に次の二つが考えられます。第一に、推定した $\lambda^ 2$ を用いてモデルの確率密度関数を計算し、それを観測された確率密度関数と比較する方法です。第二に、 $q$ を変化させたときに、 $\displaystyle{ \hat{\lambda}^ 2 _ q }$ が広い範囲でほぼ一定の値をとることを確認する方法です（詳しくは、「3. 数値実験」の部分の図を参照してください）。

■ Rで実装

　Rで実際に非ガウス性の解析をやってみます。まず、データxから $\lambda^ 2$ を推定する関数を定義します。

lambda2_est <- function(x, q = 0.25, tol = 1e-3) {
  x <- as.numeric(x)
  q <- as.numeric(q)

  # q=0, q=2 の式で現れる定数
  M_const <- (-digamma(1) + log(2)) / 2

  # 出力を事前確保
  out <- rep(NA_real_, length(q))

  # 使い回す量（速度・安定性のため）
  ax  <- abs(x)
  lax <- log(ax)

  # q≈0 判定
  is_q0 <- abs(q - 0) <= tol
  if (any(is_q0)) {
    out[is_q0] <- -mean(lax, na.rm = TRUE) - M_const
  }

  # q≈2 判定
  is_q2 <- abs(q - 2) <= tol
  if (any(is_q2)) {
    out[is_q2] <- mean(x^2 * lax, na.rm = TRUE) + M_const - 1
  }

  # 一般の q（ただし q=0,2 を除く）
  is_gen <- !(is_q0 | is_q2)
  if (any(is_gen)) {
    qg <- q[is_gen]

    # 定数部（q ごとに変わる）
    c0 <-  (log(pi) / 2)
    c1 <- -(log(2) * qg / 2)
    c2 <- -lgamma((qg + 1) / 2)

    # E[|x|^q] を q ごとに計算（q がベクトルなので sapply）
    # ※ x が長い場合ここが支配的。必要なら後述の高速化版も可能。
    m_absq <- vapply(qg, function(qq) mean(ax^qq, na.rm = TRUE), numeric(1))

    out[is_gen] <- 2 / (qg * (qg - 2)) * (c0 + log(m_absq) + c1 + c2)
  }

  # 返り値を見やすく（q の値を名前に）
  if (length(q) > 1) names(out) <- paste0("q=", format(q, digits = 6))

  # q が1つならスカラーで返す（従来互換）
  if (length(out) == 1) out[[1]] else out
}

　この関数を一度実行しておけば、以下のようにして $\displaystyle{ \hat{\lambda}^ 2 _ q }$ を計算できます。このRスクリプトでは、数値的に生成した x からパラメタを推定しています。

N <- 10000                 # サンプルサイズ
lambda2_true <- 0.4        # 真のパラメタ λ^2

# 対数正規 × 正規モデルによる擬似データの生成
x <- rnorm(N) * exp(rnorm(N, -lambda2_true, sqrt(lambda2_true)))

############################################################
# λ^2 推定の実行例
############################################################

# (1) q をスカラーで与える場合
#     単一のモーメント次数 q=0.5 を用いて λ^2 を推定
lambda2_est(x, 0.5)

# (2) q をベクトルで与える場合
#     q=0, 0.25, 0.5, 2 の各次数について λ^2 を同時に推定
#     → 戻り値は q ごとの推定値を並べたベクトル
lambda2_est(x, c(0, 0.25, 0.5, 2))

3. 数値実験

　非ガウスパラメタ $\lambda^ 2$ を正しく推定するには、ある程度長いデータが必要です。以下の図では、数値的に相乗対数正規モデルに従う乱数を生成して、確率密度関数と $\lambda^ 2$ を推定しています。

上の図はデータ数 N が10,000の場合の数値実験結果の一例です。左図に示した確率密度関数の推定結果を見ると、分布の中央（0 付近）の形状は理論分布（青破線）とよく一致しています。一方で、裾の部分についてはデータ数が十分でないため、推定の信頼性は高くありません。また、右図に示した $\displaystyle{ \hat{\lambda}^ 2 _ q }$ の $q$ 依存性については、 $\displaystyle{ q \lt 2 }$ の領域では推定値は真の値（青破線）に近いものの、 $\displaystyle{ q \gt 2 }$ の領域では大きくずれる場合があることが分かります。
　したがって、数万点以下のデータから $\lambda^ 2$ を推定する場合は、 $q$ を 0 に近い値に設定したほうが、より安定で精度の高い推定が可能です。実際には、私は $q=0.25$ 程度、あるいは $q=0$ を用いることが多いです。

4. まとめ

　今回は相乗対数正規モデルを紹介しました。このモデルの真価は、間欠性あるいは、分散の非一様性と呼ばれる現象を特徴づけるときに発揮されます。次回は、そのような現象を特徴づけるための多重スケール確率分布解析を紹介します。専門的に研究する人は以下の論文を参考にしてください。

https://journals.aps.org/pre/abstract/10.1103/PhysRevE.76.041113

https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.5.043157

https://www.sciencedirect.com/science/article/pii/016727899090035N

5. 付録：確率密度関数を描くRスクリプト

乗法log-normalモデルの確率密度関数の描画

############################################################
# 乗法log-normalモデルによる
# 非ガウス分布（確率密度関数）の数値計算と可視化
#
# ・log-normal に分布する分散 σ^2 を重ね合わせたモデル
# ・非ガウス性はパラメタ λ^2 で制御
# ・ガウス分布との比較を線形軸／対数軸で表示
############################################################


# ==========================================================
# 1) 非ガウスパラメタの設定
# ==========================================================
# λ^2 : 非ガウス性の強さを表すパラメタ
# λ^2 = 0 でガウス分布に一致し，
# λ^2 が大きいほど裾の重い分布になる
lambda.sq <- 0.4
lambda    <- sqrt(lambda.sq)


# ==========================================================
# 2) 被積分関数の定義
# ==========================================================
# 乗法ログノーマルモデルの確率密度：
#   P(x) = ∫ dσ  p_LN(σ)  N(x | 0, σ^2)
#
# ・p_LN(σ) : 分散 σ^2 が従うログノーマル分布
# ・N(x|0,σ^2) : 平均 0，分散 σ^2 の正規分布
#
# ※ xi は外側のループで与えられる評価点
multiLogNorm <- function(sig) {
  1 / (2 * pi) / lambda / sig^2 *
    exp(-(log(sig) + lambda^2)^2 / (2 * lambda^2)) *
    exp(-xi^2 / (2 * sig^2))
}


# ==========================================================
# 3) x の範囲と数値積分の設定
# ==========================================================
# 分布を評価する x の最大値
x.max <- 10

# 数値積分の検証用（ここでは直接は使用しない）
N <- 100000


# ==========================================================
# 4) 確率密度関数の数値計算
# ==========================================================
# x ≥ 0 のみ計算し，左右対称性を用いて拡張
x.list <- seq(0, x.max, length.out = 200)
n.x    <- length(x.list)

# PDF の格納ベクトルを事前確保
pdf <- numeric(n.x)

# 各 x に対して σ について数値積分
for (i in seq_len(n.x)) {
  xi     <- x.list[i]
  pdf[i] <- integrate(multiLogNorm, lower = 0, upper = Inf)$value
}

# 左右対称な分布として結合
x   <- c(-rev(x.list), x.list)
pdf <- c(rev(pdf), pdf)

# 描画用の最小値・最大値
y.min <- min(pdf)
y.max <- max(pdf)


# ==========================================================
# 5) プロット用タイトル（数式表記）
# ==========================================================
# main に λ^2 の値を数式として表示
main.txt <- bquote(lambda^2 == .(lambda.sq))


# ==========================================================
# 6) 描画（線形スケール／対数スケール）
# ==========================================================
par(mfrow = c(1, 2), mar = c(5, 5, 3, 2))


# ----------------------------------------------------------
# (a) 線形スケール
# ----------------------------------------------------------
plot(
  x, pdf, type = "l",
  col = 2, lwd = 3,
  xlim = c(-x.max, x.max),
  ylim = c(0, y.max * 1.1),
  xlab = "x", ylab = "Density",
  cex.lab = 1.7,
  main = main.txt
)

# 比較用：標準正規分布
curve(dnorm(x), add = TRUE,
      col = gray(0.5), lty = 2, lwd = 2)

legend(
  "topright",
  legend = c("Log-normal model", "Gaussian"),
  col    = c(2, gray(0.5)),
  lwd    = c(3, 2),
  lty    = c(1, 2),
  bty    = "n",
  cex    = 1.2
)


# ----------------------------------------------------------
# (b) 対数スケール（裾の違いを強調）
# ----------------------------------------------------------
plot(
  x, pdf, type = "l", log = "y",
  col = 2, lwd = 3,
  xlim = c(-x.max, x.max),
  ylim = c(y.min, y.max * 2),
  xaxt = "n", yaxt = "n",
  xlab = "x", ylab = "Density",
  cex.lab = 1.7,
  main = main.txt
)

curve(dnorm(x), add = TRUE,
      col = gray(0.5), lty = 2, lwd = 2)

# x 軸
axis(1, cex.axis = 1.2)

# y 軸（対数目盛：10^n 表記）
axis(2, at = 10^(-5:5) %x% (1:9),
     label = FALSE, tck = -0.01)

a123  <- paste(collapse = ",",
               paste0("expression(10^", -5:5, ")"))
a123b <- paste0(
  "axis(2, las=1, at=10^(-5:5), label=c(",
  a123,
  "), cex.axis=1.2, tck=-0.02)"
)
eval(parse(text = a123b))

legend(
  "topright",
  legend = c("Log-normal model", "Gaussian"),
  col    = c(2, gray(0.5)),
  lwd    = c(3, 2),
  lty    = c(1, 2),
  bty    = "n",
  cex    = 1.2
)

※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。