ケィオスの時系列解析メモランダム

時系列解析、生体情報解析などをやわらかく語ります

【Rで非ガウス統計1】正規性の検定と可視化:正規分布からのずれを読み解く

統計検定でよく耳にする 「正規性」 とは、データの分布の形が「正規分布(ガウス分布)」に従う(近似できる)という意味です。多くの統計検定や統計モデルでは、正規性が前提条件として組み込まれています。これは、正規分布が単に「よく現れる分布」だからではなく、統計理論上きわめて扱いやすい性質をもっているためです。
 もしデータ(あるいは残差)が正規分布から大きく外れているにもかかわらず、これらの手法をそのまま適用すると、

  • p 値が過小または過大に評価される
  • 有意差があるように見えて実は偶然の変動にすぎない
  • 推定された信頼区間の意味が曖昧になる

といった問題が生じる可能性があります。こうした誤用を防ぐために、解析に先立って正規性が成り立っているかどうかを確認する手段として、「正規性の検定」があります。正規性検定は、

「このデータは、統計手法が想定している前提条件を満たしているか」

を客観的にチェックするための道具です。
 もっとも、正規性検定は万能ではありません。サンプル数が大きい場合には、実務上ほとんど問題にならない程度のわずかな歪みであっても「正規分布ではない」と判定されることがあります。一方で、サンプル数が小さい場合には、本来は非正規なデータであっても、正規性があると誤って判定されることもあります。
 とくに注意が必要なのが 大規模データの場合 です。サンプルサイズが非常に大きくなると、正規分布からのごくわずかなずれであっても、検定では統計的に「正規分布である確率は非常に低い」と判定されてしまいます。これは検定の性質によるもので、データ数が増えるほど検出力(power)が高まり、実務上はほとんど無視できる程度の歪みや裾の違いまで検出されてしまうためです。
 その結果、p 値だけを見て機械的に判断すると、正規性が実質的には十分成り立っているにもかかわらず、それを棄却してしまうといった事態が起こり得ます。
 このような背景から、正規性検定の結果をそのまま信じるのではなく、「どこが、どのようにずれているのか」 を理解することが重要になります。分布の中心付近がずれているのか、裾だけが重いのか、あるいは歪度が存在するのかによって、解析結果への影響は大きく異なるからです。そのため実務では、正規性検定はあくまで参考情報と位置づけ、QQ プロットや worm プロットといった 視覚的確認法を併用して、ずれの構造を把握する という姿勢が重要になります。
 そこで本記事では、統計解析の前提条件を確認するための基本として、代表的な正規性検定の考え方と実行例、および 正規性を視覚的に確認する方法 を紹介します。検定結果の数値だけに頼るのではなく、分布の形を実際に見て「どこが、どの程度ずれているのか」を理解することを重視し、QQ プロットや worm プロットといった可視化手法についても解説します。
 ここで私が本当に伝えたいメッセージは、広く知られている一般的な方法だけでは、正規分布からのずれを定量的に評価し、その本質を理解することは難しいという点です。とくにデータ数が大規模な場合には、正規分布からのごくわずかなずれの中に、重要な情報が潜んでいることを忘れてはいけません。そのような微妙なずれは、一般的な正規性検定では捉えきれないことが多く、また QQ プロットや worm プロットといった代表的な可視化手法を用いても、必ずしも明確に見極められるとは限りません。既存の方法では見えないものがあるのです。だからこそ、私たちは既存の手法に頼るだけでなく、データの特性や研究目的に即した新しい方法を自分たちで考え、開発していく必要があります。

1. 代表的な正規性検定と R での実行例

 正規性検定にはいくつかの種類がありますが、いずれも「データが正規分布からどの程度ずれているか」を数値として評価することを目的としています。ここでは、R で広く使われている代表的な正規性検定を取り上げ、その考え方と実行例を紹介します。
 以下の検定では、

  • 帰無仮説 H _ 0:データは正規分布に従う
  • 対立仮説 H _ 1:正規分布に従わない

とし、p 値がある程度小さい場合には、「正規分布に従う」という仮定を棄却します。

1.1 Shapiro–Wilk 検定(最も標準的な正規性検定)

 Shapiro–Wilk 検定は、現在もっとも標準的に用いられている正規性検定の一つです。特に 小〜中規模のサンプルサイズ に対して高い検出力を持つことが知られており、おおよそ、 N\approx 10 から数千(R では 5000 以下)程度 のデータに対して広く用いられています。R にも標準で実装されており、多くの用途に手軽に使える正規性の検定です。

■ 基本的な考え方

 Shapiro–Wilk 検定の中心的なアイデアは非常に直感的です。「データを小さい順に並べたとき、その並び方が正規分布の並び方とどれだけ似ているか」を評価します。
 正規分布から無作為にデータを取り出した場合、その順序統計量(小さい値から順に並べたもの)は、理論的にほぼ決まった形になります。Shapiro–Wilk検定では、

  • 観測されたデータの順序統計量
  • 正規分布を仮定したときの理論的な順序統計量

一致度(相関のような量) を指標として用います。この一致度を表す統計量が W で、

  • W が 1 に近いと、正規分布に非常によく一致
  • W が小さいと、正規分布からのずれが大きい

と解釈されます。

■ R での実行例
x <- rnorm(100)

shapiro.test(x)

出力例(要旨):

Shapiro-Wilk normality test
W = 0.99, p-value = 0.73

この例では p 値が大きいため、正規性は棄却されません。すなわち、「正規分布からのずれが統計的に検出されるほど大きくはない」と判断されます。
 ただし、Shapiro–Wilk 検定を解釈する際には注意が必要です。サンプル数が非常に大きい場合には、正規分布からのごくわずかなずれであっても検定では有意と判定されやすくなります。一方で、サンプル数が極端に小さい場合には検出力が十分でなく、本来は非正規なデータであっても、そのずれを検出できないことがあります。そのため、Shapiro–Wilk 検定は正規性を確認するための 最初のチェックとしては有用 ですが、検定結果だけに依存するのではなく、必ず QQ プロットなどの視覚的確認法と併用して判断することが推奨されます。

1.2 Lilliefors 検定(Kolmogorov–Smirnov 検定の改良)

 Lilliefors 検定は、Kolmogorov–Smirnov(K–S)検定を 正規性検定として実用的に使えるよう改良した方法 です。K–S 検定では理論分布のパラメータ(平均・分散)が既知であることを前提としますが、Lilliefors 検定では 平均と分散をデータから推定した上で 正規性を評価できる点が特徴です。
 そのため、実データ解析においては、K–S 検定そのものよりも Lilliefors 検定の方が現実的な選択肢 として用いられます。

■ 基本的な考え方

 Lilliefors 検定の基本的な考え方は K–S 検定と同じで、

「観測データから得られる経験分布関数と、正規分布の累積分布関数との差を評価する」

というものです。
 手順としては、

  1. データから平均と標準偏差を推定する
  2. その推定値を用いて正規分布の累積分布関数を作る
  3. 経験分布関数との 最大差 を統計量として評価する

という流れになります。この最大差が大きいほど、分布のどこかに「正規分布からのずれ」が存在すると判断されます。
 ただし、平均・分散をデータから推定しているため、通常の K–S 検定の分布は使えず、Lilliefors 検定専用に導出された近似分布やシミュレーション結果に基づいて p 値が計算されます。

■ R での実行例

 Lilliefors 検定は base R には含まれていないため、nortest パッケージを用います。

install.packages("nortest")
library(nortest)

x <- rnorm(100)

lillie.test(x)

出力例(要旨):

Lilliefors (Kolmogorov-Smirnov) normality test
D = 0.05, p-value = 0.62

この例では p 値が大きいため、正規性は棄却されません。

■ 実務での位置づけと注意点

 Lilliefors 検定は、K–S 検定の「パラメータ既知」という非現実的な前提を取り除いた点で、正規性検定としては理論的に妥当な改良と言えます。一方で、

  • 分布の中心付近のずれに主に反応する
  • 裾の違いには必ずしも敏感ではない
  • サンプルサイズが大きいと、やはりわずかなずれでも有意になりやすい

といった性質は残っています。そのため、Lilliefors 検定も 単独で正規性を判断するための決定打ではなく、 Shapiro–Wilk 検定や QQ プロットなどと併用し、「どこがどのようにずれているのか」を確認するための補助的な手段として使うのが適切です。

1.3 Anderson–Darling 検定(裾を重視)

 Anderson–Darling(A–D)検定は、分布の裾(tail)のずれを重視して正規性を評価する検定です。正規分布からのわずかなずれであっても、とくに極端な値(外れ値や裾の重さ)が重要になる解析において、有効な正規性検定として用いられます。Shapiro–Wilk検定や Lilliefors検定が分布全体の一致度を比較的均等に評価するのに対し、Anderson–Darling検定は 分布の端の部分により強い重みを与える という点が大きな特徴です。

■ 基本的な考え方

 Anderson–Darling 検定も、基本的な枠組みは

「観測データの分布が、理論的な正規分布とどれだけ一致しているか」

を評価するものです。ただし、累積分布関数(CDF)を用いる点は Lilliefors 検定と同じでも、Anderson–Darling 検定では、

  • 分布の中心付近よりも
  • 裾(非常に小さい値・非常に大きい値)に大きな重み

を与えてずれを評価します。その結果、

  • 裾が重い分布
  • 外れ値を含む分布
  • 極端値の頻度が正規分布と異なる場合

に対して、ずれを検出しやすくなっています。この重み付けによって計算される統計量が A^ 2(エースクエア)で、

  • A^ 2 が小さいと、正規分布によく一致
  • A^ 2 が大きいと、裾を含めて正規分布からのずれが大きい

と解釈されます。

■ R での実行例

 Anderson–Darling 検定も、base R には含まれていないため、nortest パッケージを用います。

install.packages("nortest")
library(nortest)

x <- rnorm(100)

ad.test(x)

出力例(要旨):

Anderson-Darling normality test
A = 0.29, p-value = 0.61

この例では p 値が大きいため、正規性は棄却されません。

■ 実務での位置づけと注意点

 Anderson–Darling 検定は、裾の違いが解析結果に大きく影響する場合 に特に有用です。たとえば、

  • 外れ値の影響が問題になる解析
  • リスク評価や極端現象の解析
  • 正規分布の裾の軽重が重要な意味を持つ場合

などでは、Shapiro–Wilk検定よりも有用な情報を与えることがあります。
 一方で、Anderson–Darling 検定もまた、

  • サンプルサイズが大きいと、わずかな裾の違いでも有意になりやすい
  • どの部分がどのようにずれているかまでは直接示さない

という性質を持っています。そのため、Anderson–Darling検定も 単独で結論を出すための検定ではなく、 QQプロットや wormプロットなどの視覚的確認法と併用し、「裾のどこに、どのようなずれがあるのか」を理解するための補助的な手段として使うのが適切です。

2. 正規性の視覚的確認法

 前節では、Shapiro–Wilk 検定、Lilliefors 検定、Anderson–Darling 検定といった、代表的な正規性検定を紹介しました。これらの検定は、正規分布からのずれを 数値として客観的に評価できる という点で有用ですが、検定結果だけで正規性の妥当性を判断することには限界があります。
 とくに、サンプルサイズが大きい場合には、実務上ほとんど問題にならない程度のわずかなずれであっても「有意」と判定されることがあります。一方で、サンプルサイズが小さい場合には、明らかな非正規性があっても検出できないこともあります。このように、正規性検定の結果は データ数に強く依存する ため、p 値だけを基準に機械的な判断を行うことは適切ではありません。
 そこで重要になるのが、分布の形を直接「目で確認する」視覚的確認法 です。視覚的な方法を用いることで、

  • 正規分布からのずれがどの部分に現れているのか
  • 中心付近のずれなのか、裾の違いなのか
  • 歪度や外れ値がどの程度存在するのか

といった情報を直感的に把握することができます。これは、検定統計量や p 値からは直接読み取れない重要な情報です。
 以下では、正規性の視覚的確認法として広く用いられている QQ プロット と、QQ プロットのずれをより明確に可視化する worm プロット を取り上げ、それぞれの考え方と R による実行例を解説します。正規性検定の結果とこれらの可視化を組み合わせることで、正規分布の仮定が解析目的に対して妥当かどうかを、より適切に判断できるようになります。

2.1 QQ プロット(Quantile–Quantile プロット)

 QQ プロット(Quantile–Quantile プロット)は、データの分布が正規分布とどの程度一致しているかを視覚的に確認するための最も基本的な方法です。正規性検定が「ずれの大きさ」を数値として評価するのに対し、QQ プロットは 「どこが、どのようにずれているのか」 を直感的に示してくれます。

■ 基本的な考え方

QQ プロットでは、

  • 横軸:正規分布を仮定したときの理論分位点(theoretical quantiles)
  • 縦軸:観測データの分位点(sample quantiles)

を対応させてプロットします。正規分布に従うデータであれば、観測分位点と理論分位点はほぼ 1 対 1 に対応するため、点は一直線上に並びます。この直線からのずれが、そのまま「正規分布からのずれ」を意味します。
 Shapiro–Wilk 検定が「並び方の一致度」を数値化しているのに対し、QQ プロットは その並び方を可視化したもの と考えると理解しやすいでしょう。

■ R での実行例

 R では、QQ プロットは非常に簡単に描くことができます。

x <- rnorm(100)

qqnorm(x)
qqline(x, col = "red")

qqnorm() で QQ プロットを描き、qqline() で基準となる直線を重ねています。

■ QQ プロットの読み方

 QQ プロットでは、点の並び方から正規分布からのずれの特徴を読み取ります。

  • 点がほぼ直線上に並んでいれば、正規分布とよく一致している

  • S 字状に曲がると、分布に歪度(skewness)がある

  • 両端が大きく外側に膨らむと、裾が重い(heavy tail)

  • 両端が内側に入り込むと、裾が軽い(light tail)

このように、QQ プロットは、「正規かどうか」だけでなく、「どのように正規から外れているか」を把握するために非常に有用です。

■ 正規性検定との関係

 正規性検定では p 値という一つの数値しか得られませんが、QQ プロットでは、

  • ずれが中心付近に集中しているのか
  • 裾だけに現れているのか
  • 外れ値が原因なのか

といった情報を直接確認できます。そのため実務では、

  1. QQ プロットで分布の形を確認する
  2. 正規性検定の結果を参考情報として解釈する

という順序で判断することが推奨されます。

次節では、QQ プロットのずれをさらに明確にするための手法として、worm プロット を紹介します。

2.2 worm プロット

 worm プロットは、QQ プロットから基準となる直線を差し引いた図であり、QQ プロットでは見えにくい 微妙なずれの構造を強調して可視化する方法です。その形が「うねるミミズ(worm)」のように見えることから、この名前が付けられています。
 QQ プロットが「分布全体の一致度」を直感的に示すのに対し、worm プロットは 正規分布からの偏差そのものを直接描く 点に特徴があります。

■ 基本的な考え方

 worm プロットの構造は単純で、

  • 横軸:正規分布を仮定したときの理論分位点
  • 縦軸:観測分位点 − 理論分位点

をプロットします。すなわち、QQ プロット上での「直線からのずれ」だけを取り出して描いたもの が worm プロットです。正規分布に従うデータであれば、理論分位点と観測分位点の差はほぼ 0 となり、点は 横軸(0 の直線)付近にランダムに散らばります
 一方で、正規分布からの系統的なずれがある場合には、その構造が明確なパターンとして現れます。

正規乱数のQQプロット(左)と wormプロット(右)。QQプロットでは、点が赤線上にあること、wormプロットでは、中央で水平になり、上下の破線を超えないことを確認する。

■ R での実行例

 以下のRスクリプトで、wormプロットを描けます。この例では、QQプロットも同時に描いています。

# ----------------------------------------------------------
# データ生成(自分のデータがある場合は、コメントアウトしてください)
# ----------------------------------------------------------
x <- rnorm(100)
n <- length(x)

# ----------------------------------------------------------
# 計算パート(Worm plot用)
# ----------------------------------------------------------
# 1) 前処理:標準化と並べ替え
x.stdzd <- as.numeric(scale(x))
x.stdzd.o <- sort(x.stdzd)

# 2) 理論分位点
qt.plot <- function(k, n) (k - 0.5) / n
qt <- qnorm(p = qt.plot(1:n, n))

# 3) 偏差(Deviation)の計算
sd.x <- sd(x.stdzd)
x.detrend <- x.stdzd.o - sd.x * qt

# 4) 95% 参考帯(±1.96 SE)の計算
z.SE <- function(z, n) {
  sqrt(pnorm(z) * (1 - pnorm(z)) / n) / dnorm(z)
}

# 描画範囲設定
x.abs.max <- ceiling(max(abs(x.stdzd), na.rm = TRUE) / 0.5) * 0.5
x.range   <- c(-x.abs.max, x.abs.max)

# 滑らかな参考帯用のグリッド
z.plot <- seq(-x.abs.max, x.abs.max, length.out = 500)
SE <- sd.x * z.SE(z.plot, n)

# ----------------------------------------------------------
# プロット設定
# ----------------------------------------------------------
# 画面分割とマージン設定(余白を少し調整してタイトルやラベルが見切れにくくする)
par(mfrow = c(1, 2), mar = c(5, 5, 3, 2), family = "sans")

# 共通のデザインパラメータ
pt.col   <- rgb(0, 0.4, 0.8, 0.7)  # 少し透明な青
line.col <- "firebrick"            # 落ち着いた赤
ci.col   <- "gray20"               # 信頼区間用
pt.cex   <- 1.2                    # 点のサイズ
lab.cex  <- 1.3                    # ラベルサイズ
axis.cex <- 1.1                    # 軸の数字サイズ

# ----------------------------------------------------------
# 左:Normal Q-Q Plot
# ----------------------------------------------------------
# プロット(データを描画せず枠だけ作る)
qq <- qqnorm(x, plot.it = FALSE)
plot(qq$x, qq$y, type = "n",
     main = "Normal Q-Q Plot",
     xlab = "Theoretical Quantiles",
     ylab = "Sample Quantiles",
     cex.lab = lab.cex, cex.axis = axis.cex, las = 1)

# グリッドを追加
grid(lwd = 1, col = "gray90")

# 点と線を描画
points(qq$x, qq$y, pch = 16, col = pt.col, cex = pt.cex)
qqline(x, col = line.col, lwd = 2)

# 凡例 (左上)
legend("topleft", 
       legend = c("Observed", "Theoretical"),
       col = c(pt.col, line.col),
       pch = c(16, NA), 
       lty = c(NA, 1), 
       lwd = c(NA, 2),
       bty = "n", cex = 1)

# ----------------------------------------------------------
# 右:Worm Plot (Detrended Q-Q Plot)
# ----------------------------------------------------------
# Y軸の範囲を動的に計算(CIかデータの大きい方に合わせる)
max.y <- max(c(abs(x.detrend), 1.96 * SE), na.rm = TRUE)
ylim.dynamic <- c(-max.y, max.y) * 1.1 # 少し余裕を持たせる

plot(qt, x.detrend,
     type = "n", # まだ描画しない
     xlab = "Unit Normal Quantile",
     ylab = "",
     xlim = x.range, 
     ylim = c(-1.2, 1.2),
     xaxs = "i", yaxs = "i",
     main = "Worm Plot",
     cex.lab = lab.cex, cex.axis = axis.cex, las = 1)

# グリッド
grid(lwd = 1, col = "gray90")

# Y軸ラベル(回転させる)
mtext("Deviation", side = 2, line = 3.5, cex = lab.cex, las = 0)

# 信頼区間 (CI) の描画
lines(z.plot,  1.96 * SE, col = ci.col, lwd = 1.5, lty = 2)
lines(z.plot, -1.96 * SE, col = ci.col, lwd = 1.5, lty = 2)

# 基準線 (0線)
abline(h = 0, col = "gray50", lwd = 1, lty = 1)
abline(v = 0, col = "gray50", lwd = 1, lty = 3) # 中心線は見やすく点線で

# データの点
points(qt, x.detrend, pch = 16, col = pt.col, cex = pt.cex)

# 凡例 (右上または右下、データと被らない場所が良いがここでは右上に固定)
legend("topright", 
       legend = c("Deviation", "95% CI Limit"),
       col = c(pt.col, ci.col),
       pch = c(16, NA), 
       lty = c(NA, 2), 
       lwd = c(NA, 1.5),
       bty = "n", cex = 1)

ここでは、

  • theo:正規分布の理論分位点
  • obs:観測データの分位点

を計算し、その差をプロットしています。

■ worm プロットの読み方

 worm プロットでは、偏差の形から正規分布からのずれの特徴を読み取ります。

  • 点が 0 の周りにランダムに散らばるのであれば、正規分布とよく一致している

  • 中央付近が上または下に曲がると、歪度(skewness)が存在する

  • 両端が大きく正または負に振れると、裾が重い、または外れ値の影響がある

  • 系統的な曲線構造が現れると、正規分布からの構造的なずれが存在する

QQ プロットでは「ほぼ直線」に見えても、worm プロットでは 小さなずれが明確なパターンとして現れる ことがあります。

■ QQ プロットとの使い分け

 QQ プロットは分布の全体像を把握するのに適しており、 worm プロットは そのずれを詳細に分析するための拡大鏡 のような役割を果たします。  実務では、

  1. QQ プロットで大まかな分布形状を確認する
  2. worm プロットでずれの構造を詳しく調べる

という順序で用いると、正規性の判断がより明確になります。

3. 非正規性を示すデータ生成

 正規性検定や QQ プロットを理解するうえで重要なのは、「非正規」と一言で言っても、その中身は一様ではないという点です。分布が正規分布から外れる理由には、いくつか典型的なパターンがあります。
 ここでは、非正規性の代表的な 3 つのタイプを意図的に作り出すデータ生成モデルを紹介します。

  1. 歪度(非対称性)をもつ分布
  2. 裾が重くなる分布(分散ゆらぎ)
  3. 裾のみに外れ値が入る分布

これらは、正規性検定や視覚的確認法が「どの特徴に敏感なのか」を理解するための、非常に良い対照例になります。

3.1 skew-normal 分布:歪度のみをもつ非正規分布

 skew-normal 分布は、正規分布を拡張し、分布の歪度(左右非対称性)だけを連続的に制御できるモデルです。

  • 正規分布:左右対称
  • skew-normal:片側に尾が長い

という違いがありますが、裾の重さ自体は正規分布とほぼ同じです。この分布は、以下のパラメタで形を指定します。

  • \xi:位置パラメータ(中心)
  • \omega:スケールパラメータ(分散)
  • \alpha:歪みパラメータ

ここでは、

  • \alpha = 0:正規分布
  • \alpha > 0:右に歪んだ分布
  • \alpha \lt 0:左に歪んだ分布

となり、\alpha の値で非正規性を調節できます。
 以下のRスクリプトで、\alpha の値を0から、少しずつ変えて乱数xを生成し、それを使って正規性の検定やプロットをしてみてください。

############################################################
# Skew-normal Distribution
#
# モデル:
#   X ~ SN(xi, omega, alpha)
#
#   alpha  : 歪みパラメータ
#
# 解釈:
#   - alpha = 0 → 正規分布 N(xi, omega^2)
#   - alpha > 0 → 右に歪んだ分布
#   - alpha < 0 → 左に歪んだ分布
############################################################


# ----------------------------------------------------------
# 必要パッケージの自動チェック&インストール
# ----------------------------------------------------------
# skew-normal 分布は base R には含まれていないため,
# sn パッケージを利用する。
# 未インストールの場合のみ自動的にインストール。
if (!requireNamespace("sn", quietly = TRUE)) {
  install.packages("sn")
}

library(sn)


# ----------------------------------------------------------
# skew-normal 乱数生成関数
# ----------------------------------------------------------
# n     : サンプル数
# alpha : 歪みパラメータ(0で正規分布)
# xi    : 位置パラメータ
# omega : スケールパラメータ
#
# 戻り値:
#   数値ベクトル x(skew-normal 分布に従う)
# ----------------------------------------------------------
rskewnorm <- function(n, alpha = 0, xi = 0, omega = 1, seed = NULL) {
  if (!is.null(seed)) set.seed(seed)
  
  # sn::rsn は skew-normal 分布の乱数生成関数
  x <- rsn(n, xi = xi, omega = omega, alpha = alpha)
  
  # 念のため numeric に変換
  return(as.numeric(x))
}


# ----------------------------------------------------------
# 使用例:パラメータ設定
# ----------------------------------------------------------
n     <- 1000   # データ数
alpha <- 2      # 歪みパラメータ(正 → 右に歪む)

# 乱数生成
x <- rskewnorm(n = n, alpha = alpha, seed = 1)


# ----------------------------------------------------------
# 可視化:ヒストグラム + 理論分布
# ----------------------------------------------------------
hist(x, breaks = 40, probability = TRUE,
     main = paste("skew-normal (alpha =", alpha, ")"),
     xlab = "x")

# 理論的 skew-normal の確率密度関数
curve(dsn(x, alpha = alpha),
      add = TRUE, col = "red", lwd = 2)

# 比較用:同じ平均・分散の正規分布
# (歪度がない場合との違いを明示)
curve(dnorm(x),
      add = TRUE, col = "blue", lwd = 2, lty = 2)


# ----------------------------------------------------------
# 凡例
# ----------------------------------------------------------
legend("topright",
       legend = c("Histogram (simulation)",
                  paste("Skew-normal (alpha =", alpha, ")"),
                  "Normal distribution"),
       col    = c("gray", "red", "blue"),
       lwd    = c(NA, 2, 2),
       lty    = c(NA, 1, 2),
       pch    = c(15, NA, NA),
       pt.cex = 1.5,
       bty    = "n")

QQプロット(左)とworm プロット(右)。

3.2 lognormal scale mixture:分散ゆらぎによる厚い裾分布

 次に紹介するのは、分散がランダムにゆらぐことによって裾が重くなるタイプの非正規分布です。このモデルは、lognormal scale mixture of Gaussians、multiplicative log-normal model (私がつけた名前)と呼ばれます。理屈を簡単に説明すれば、正規分布に従う標準正規乱数 Z \sim \mathcal N(0,1) と対数正規分布に従う乱数  Y = \exp(U)(ここで、U \sim \mathcal N(-\lambda^ 2,\lambda^ 2))をかけているだけです。
 このモデルでは、条件付きでは正規分布に従いますが、その分散 Y^ 2 が対数正規分布に従って変動します。その結果、周辺分布としては裾の重い分布が得られます。

パラメータ  \lambda の意味は次のとおりです。

  •  \lambda = 0:分散のゆらぎがなく、正規分布に一致します
  •  \lambda > 0:分散のゆらぎが大きくなり、裾が重くなります

 以下のRスクリプトで、\lambda の値を0から少しずつ増加させて乱数xを生成し、それを使って正規性の検定やプロットをしてみてください。

############################################################
# Lognormal Scale Mixture of Gaussians
# (multiplicative log-normal model)
# モデル:
#   Z ~ N(0, 1)
#   U ~ N(-lambda^2, lambda^2)
#   Y = exp(U)        (対数正規分布)
#   X = Y * Z
#
#  lambda : 非正規パラメタ
#
# 特徴:
#   - lambda = 0 → 純粋な正規分布
#   - lambda > 0 → 中央はほぼ正規、裾が重くなる
#   - 正規性検定では裾の影響で棄却されやすい
############################################################


# ----------------------------------------------------------
# 乱数生成関数:標準正規 × 対数正規
# ----------------------------------------------------------
r_lognorm_scaled_normal <- function(n, lambda, seed = NULL) {
  if (!is.null(seed)) set.seed(seed)
  
  # 1) 標準正規乱数 Z ~ N(0, 1)
  z <- rnorm(n)
  
  # 2) 対数正規乱数 Y = exp(U),
  #    U ~ N(-lambda^2, lambda^2)
  #    → 分散の平均的スケールを抑えた設定
  y <- rlnorm(n, meanlog = -lambda^2, sdlog = lambda)
  
  # 3) 乗算(分散がランダムに変動する正規分布)
  x <- y * z
  
  return(x)
}


# ----------------------------------------------------------
# 確率密度関数 f_X(x) を数値積分で評価
# ----------------------------------------------------------
# 理論式:
#   f_X(x) = ∫ N(x | 0, y^2) f_Y(y) dy
#          = ∫ dnorm(x, sd=y) dlnorm(y) dy
#
# ※ 閉形式は無いので数値積分を用いる
# ----------------------------------------------------------
d_lognorm_scaled_normal <- function(x, lambda) {
  
  # 対数正規分布 f_Y(y)
  fy <- function(y) {
    dlnorm(y, meanlog = -lambda^2, sdlog = lambda)
  }
  
  # 各 x について積分を実行
  fx <- sapply(x, function(xx) {
    
    # 被積分関数:条件付き正規 × mixing density
    integrand <- function(y) {
      dnorm(xx, mean = 0, sd = y) * fy(y)
    }
    
    integrate(integrand,
              lower = 0,
              upper = Inf,
              rel.tol = 1e-6)$value
  })
  
  return(fx)
}


# ----------------------------------------------------------
# 使用例:パラメータ設定
# ----------------------------------------------------------
n      <- 1000   # サンプル数
lambda <- 0.3    # 分散ゆらぎの強さ(0で正規分布)

# 乱数生成
x <- r_lognorm_scaled_normal(n, lambda, seed = 1)

# x 軸用グリッド(pdf 描画用)
xg <- seq(min(x), max(x), length.out = 400)

# 数値積分による理論 pdf
fx <- d_lognorm_scaled_normal(xg, lambda)


# ----------------------------------------------------------
# 可視化:ヒストグラム + 理論 pdf
# ----------------------------------------------------------
hist(x, breaks = 120, probability = TRUE,
     main = bquote(
       X == Y %.% Z ~
       (log~Y %~% N(-lambda^2, lambda^2)) ~ "," ~
       lambda == .(lambda)
     ),
     xlab = "x")

# 数値積分による pdf(理論)
lines(xg, fx, col = "red", lwd = 2)

# 比較用:標準正規分布
curve(dnorm(x), add = TRUE,
      col = "blue", lwd = 2, lty = 2)


# ----------------------------------------------------------
# 凡例
# ----------------------------------------------------------
legend("topright",
       legend = c("Histogram (simulation)",
                  "Numerical PDF (scale mixture)",
                  "Standard normal"),
       col    = c("gray", "red", "blue"),
       lwd    = c(NA, 2, 2),
       lty    = c(NA, 1, 2),
       pch    = c(15, NA, NA),
       pt.cex = 1.5,
       bty    = "n")

QQプロット(左)とworm プロット(右)。

3.3 裾のみ Weibull 外れ値をもつモデル

 最後に紹介するのは 「裾のみ Weibull 外れ値をもつモデル」です。これは、私が作ったモデルです。まだ、論文に仕上げていませんが、実世界の分布に意外と合うと思います。このモデルでは、分布の中心部分は正規分布に厳密に従う一方で、外れ値はあらかじめ定めた境界  |X|>T を超えた 裾の領域にのみ 出現するように設計されています。さらに、その外れ値の大きさ分布として Weibull 分布を仮定することで、裾の減衰の速さや極端値の出現頻度を柔軟に制御できる構造を持たせています。

\displaystyle{
X =
\begin{cases}
\mathcal N(\mu,\sigma^2), & \text{with probability } 1-\varepsilon, \\
\mu + \sigma \, S \,(T + W), & \text{with probability } \varepsilon,
\end{cases}
}

ここで、\displaystyle{
S=\pm1
} は外れ値が左右対称に現れることを表す符号変数であり、\displaystyle{
W\sim\mathrm{Weibull}
} は外れ値の大きさを決める確率変数です。

各パラメータの意味は次のとおりです。

  • \displaystyle{\varepsilon}:外れ値が混入する割合を表します
  • T:外れ値が出現する境界であり、必ず |X|>T となります
  • Weibull 分布の shape:裾の減衰の速さ、すなわち極端値の出やすさを制御します

 以下のRスクリプトで、\varepsilon の値を0から少しずつ増加させて乱数xを生成し、それを使って正規性の検定やプロットをしてみてください。

############################################################
# Tail-contaminated Normal Model
#
# モデル:
#   X ~ N(mu, sigma^2)                     with prob. 1 - eps
#   X = mu + sigma * S * (T + W)           with prob. eps
#     S ∈ {−1, +1}  (左右対称)
#     W ~ Weibull(shape, scale)
#
#  eps   : 外れ値の割合
#
############################################################

# ----------------------------------------------------------
# 関数:裾のみ Weibull 外れ値を持つ正規乱数生成
# ----------------------------------------------------------
r_tail_weibull_outlier_normal <- function(
  n,
  eps = 0.01,        # 外れ値が入る確率
  T = 3,             # 外れ値が入る境界 |X| > T
  shape = 1.5,       # Weibull 形状母数(小さいほど重い裾)
  scale = 1,         # Weibull スケール母数
  mu = 0,            # 正規分布の平均
  sigma = 1,         # 正規分布の標準偏差
  seed = NULL
) {
  if (!is.null(seed)) set.seed(seed)
  
  # 1) 基本となる正規乱数
  x <- rnorm(n, mean = mu, sd = sigma)
  
  # 2) 外れ値になるインデックスを決定
  is_out <- runif(n) < eps
  m <- sum(is_out)
  
  # 外れ値が無い場合はそのまま返す
  if (m == 0) return(x)
  
  # 3) Weibull 分布に従う「裾の大きさ」を生成
  #    → 必ず |X| > T になるように T を加える
  w <- rweibull(m, shape = shape, scale = scale)
  mag <- T + w
  
  # 4) 左右対称にするため符号をランダムに付与
  sgn <- sample(c(-1, 1), m, replace = TRUE)
  
  # 5) 外れ値で置き換え
  #    境界は一般には mu ± sigma*T
  x[is_out] <- mu + sigma * sgn * mag
  
  return(x)
}

# ----------------------------------------------------------
# 使用例:パラメータ設定
# ----------------------------------------------------------
n     <- 10000      # サンプル数
eps   <- 0.02        # 外れ値の割合(0なら純粋な正規分布)
T     <- 2.5         # 裾の開始位置
shape <- 1.8         # Weibull shape(<2 で重い裾)
scale <- 1           # Weibull scale

# 乱数生成
x <- r_tail_weibull_outlier_normal(
  n     = n,
  eps   = eps,
  T     = T,
  shape = shape,
  scale = scale,
  seed  = 1
)

# ----------------------------------------------------------
# 可視化:ヒストグラム + 参照線
# ----------------------------------------------------------
hist(x, breaks = 100, probability = TRUE,
     main = sprintf(
       "Tail Weibull outliers: eps=%.3f, T=%.1f, shape=%.2f",
       eps, T, shape
     ),
     xlab = "x")

# 比較用:標準正規分布
curve(dnorm(x), add = TRUE, col = "blue", lwd = 2, lty = 2)

# 外れ値が入る境界 |X| = mu ± sigma*T
abline(v =  mu + sigma * T, col = "darkgreen", lwd = 2, lty = 3)
abline(v =  mu - sigma * T, col = "darkgreen", lwd = 2, lty = 3)

# ----------------------------------------------------------
# 凡例
# ----------------------------------------------------------
legend("topright",
       legend = c("Histogram (simulation)",
                  "Standard normal",
                  expression("|X| = T (outlier boundary)")),
       col    = c("gray", "blue", "darkgreen"),
       lwd    = c(NA, 2, 2),
       lty    = c(NA, 2, 3),
       pch    = c(15, NA, NA),
       pt.cex = 1.5,
       bty    = "n")

QQプロット(左)とworm プロット(右)。

4. まとめ

 非正規分布の例のうち後半の2つは、私が研究してきたものです。多くの方々は、正規か、非正規かくらいしか気にしないかもしれませんが、私は分布を読むということを大切にして研究しています。みなさんも、分布に興味をもってもらえると嬉しいです。
 今回は説明が長くなりすぎて、本当に伝えたいメッセージが書けなかったので、続きを別で書きます。

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