【Rで非ガウス統計２】非対称なSkew Normal分布：２次元正規からの切り取り

今回は、非対称性を示す分布モデルである Skew Normal 分布を紹介します。このモデルは、一つのパラメタを用いて、正規分布からずれた非対称な分布を表現できる点が特徴です。その定義式をみれば一見難しい理論に基づいているように感じますが、その導出のアイデアはとても簡単です。
　Skew Normal 分布は、Adelchi Azzalini によって 1985 年に導入されました。

Azzalini, A. A class of distributions which includes the normal ones Scandinavian Journal of Statistics, 1985

当時すでに「非対称分布」は数多く提案されていましたが、その多くは、パラメタが多すぎて数学的に扱いにくい、あるいは 正規分布との関係が不明瞭といった問題を抱えていました。
　そのような状況の中で、Azzalini は比較的簡単な仮定の下で、正規分布を非対称にゆがめるアイデアを提案しました（下図）。彼のモデルでは、パラメタは一つだけとしながら、正規分布からのずれ（非対称性）を特徴づけること、そして、分布を生む仕組みがわかりやすいということです。仕組みといっても、難解なものではなく、正規分布の自然な拡張として理解できる点が、この分布の良い点です。

Skwe Normal分布のイラスト。右側のパネルの赤破線（右上段）と赤実線（右下段）のグラフがskew normal。左上段は２次元正規分布から、ピンク色の点だけど抽出し、skew normalが実現する仕組みの説明。 $\alpha=0$ のとき正規分布になり、 $\alpha \neq 0$ で非対称性が生じる。

Skwe Normal分布のイラスト。右側のパネルの赤破線（右上段）と赤実線（右下段）のグラフがskew normal。左上段は２次元正規分布から、ピンク色の点だけど抽出し、skew normalが実現する仕組みの説明。 $\alpha=0$ のとき正規分布になり、 $\alpha \neq 0$ で非対称性が生じる。

1. Skew Normal の定義
2. 非対称性を生むアイデア：2 次元正規分布から一部を切り取る
3. RでSkew Normal
4. まとめ

1. Skew Normal の定義

　式だけを見ると、少し難解に見えるかもしれませんが、まずは定義式を書いておきます。なぜこの形になるのかは、後で説明します。
　Skew Normal 分布は次のような確率密度関数で定義されます。

$\displaystyle{ f(x)= 2\, \phi(x)\, \Phi(\alpha x) }$

ここで、

$\displaystyle{ \phi(x) = \frac{1}{\sqrt{2\pi}} \exp\!\left(-\frac{x^{2}}{2}\right) }$

は、平均 0、分散 1 の 標準正規分布の確率密度関数、

$\displaystyle{ \Phi(x) = \int_{-\infty}^{x} \phi(t)\,dt }$

は、その 分布関数（誤差関数）を表します。したがって、Skew Normal 分布の確率密度関数を積分を含む形で書けば、

$\displaystyle{ f(x) = \frac{1}{\pi} \exp\!\left(-\frac{x^{2}}{2}\right) \int_{-\infty}^{\alpha x} \exp\!\left(-\frac{t^{2}}{2}\right)\,dt }$

となります。
　このモデルで非対称性をコントロールするのは、パラメタ $\displaystyle{ \alpha \in \mathbb{R} }$ で、

$\alpha=0$ のとき、 Skew Normal 分布は正規分布に一致し、
$\displaystyle{\alpha\gt0}$ では右側の裾が厚い非対称分布、
$\displaystyle{\alpha\lt0}$ では左側の裾が厚い非対称分布、

となります。 $\alpha=0$ のとき、 $\displaystyle{ \Phi(0)= \int _ {-\infty}^ {0} \phi(t)\,dt = 1/2 }$ ですので、 $f(x)= \phi(x)$ となり標準正規分布になります。
　一見すると、「なぜこのような形になるのだろうか」と疑問に感じる方も多いと思います。しかし、Skew Normal 分布は 2 次元正規分布からの条件付き分布として理解することができ、その導出過程をたどると、数式としては納得のいく構成になっています。一方で、このようなメカニズムが自然界や実データの背後にどの程度普遍的に存在しているのかについては、正直なところ、私には判断できません。

2. 非対称性を生むアイデア：2 次元正規分布から一部を切り取る

　Skew Normal 分布を生成過程に基づいて理解するポイントは、「1 次元の分布を直接いじる」のではなく、「2 次元に拡張してから一部を切り取る」ということです。
　まず、独立な 2 つの標準正規乱数 $X$ と $Y$ を考えます。

$\displaystyle{ (X, Y) \sim N_2(\mathbf{0}, I) }$

ここで、 $\sim$ は、確率分布が同じという意味で、 $I$ は２次元単位行列です。ようするに、 $X$ と $Y$ は、互いに独立な、平均ゼロ、分散１の正規乱数ということです。 $X$ と $Y$ をxy平面に描くと、これは、下の図左のような平面上に広がる円対称なガウス雲になります。この段階では、当然ながら $X$ の周辺分布は完全に対称な正規分布です。

　　　ここで、Azzalini は次のような、非常にシンプルな操作を導入しました。

すべての点を使うのではなく、ある条件を満たした点だけを使う。

具体的には、パラメタ $\alpha$ を導入し、次の条件を考えます。

$\displaystyle{ Y \le \alpha X }$

これは、上の図右のように、直線 $y=\alpha x$ の下側の半平面を表しています。 $\alpha$ は直線の傾きで、この傾きが分布の歪みを決めます。
　ここで定義する確率変数は、

$\displaystyle{ Z = X \mid (Y \le \alpha X) }$

です。ここで、「 $\mid$ 」は、「右側に書いた条件の下で」という意味です。すなわち、「2 次元正規分布から点を生成し、そのうち条件「 $Y \le \alpha X$ 」を満たした点の $X$ 成分だけを見る」という操作を行っています。
　このとき得られる $Z$ の分布が、Skew Normal 分布です。このときの密度関数は、条件付き分布の定義から

$\displaystyle{ f_Z(x) = \frac{P(Y \le \alpha x \mid X = x)\, f_X(x)} {P(Y \le \alpha X)} }$

と書けます。ここで、 $X$ と $Y$ は 独立な標準正規分布として生成されているため、「 $X=x$ の条件」を固定しても、何の影響もないので、

$\displaystyle{ P(Y \le \alpha x \mid X = x) = P(Y \le \alpha x) }$

と条件を消すことができます。さらに、 $Y$ は正規分布、

$\displaystyle{ Y \sim N(0,1) }$

に従うので、

$\displaystyle{ P(Y \le \alpha x) = \int_{-\infty}^{\alpha x} \frac{1}{\sqrt{2 \pi}} e^{- \frac{t^2}{2}}\,dt = \Phi(\alpha x) }$

となります。これが、分布の定義式に、誤差関数 $\Phi(\alpha x)$ が現れる理由です。
　分母 $P(Y \le \alpha X)$ については、2 次元正規分布は原点に関して対称であるため（分布を均等に２つに分けるため）、

$\displaystyle{ P(Y \le \alpha X) = \frac{1}{2} }$

となります。
　以上をまとめると、Skew Normal 分布の確率密度関数として、

$\displaystyle{ f_Z(x)= 2\,\phi(x)\,\Phi(\alpha x) }$

が得られます。以上を何となく理解したうえで、最初に掲載したグラフアニメを見て、理解とイメージをつなげてください。

3. RでSkew Normal

　Skew Normal 分布を R で扱うための方法を紹介しておきます。Skew Normal を含む Azzalini 系の分布（skew-normal / skew-t など）を扱うには、sn パッケージを使うのが簡単です。

3.1 パッケージの読み込み

　sn パッケージを既にインストールしてある人は、この部分は必要ありませんが、sn が未インストールの可能性がある人は以下のスクリプトを実行して下さい。


# ============================================================
# 0) Utility: install + load packages automatically
# ============================================================
load_or_install <- function(pkgs) {
  for (p in pkgs) {
    if (!requireNamespace(p, quietly = TRUE)) {
      install.packages(p)
    }
    library(p, character.only = TRUE)
  }
}

# Use 'sn' for skew-normal / skew-t, etc.
load_or_install(c("sn"))

3.2 確率密度関数（PDF）の描画

　sn パッケージでは、Skew Normal の密度関数は dsn() で計算できます。パラメータは

xi：位置（location）
omega：尺度（scale）
alpha：歪み（shape）

です。alpha=0 のとき、Skew Normal は正規分布に一致します。


# ============================================================
# 1) PDF plot (Skew-normal vs Normal)
# ============================================================
xi    <- 0
omega <- 1
alpha <- 3

x <- seq(-4, 4, length.out = 1000)

pdf_skew <- dsn(x, xi = xi, omega = omega, alpha = alpha)

plot(x, pdf_skew, type = "l", lwd = 2,
     xlab = "x", ylab = "density",
     main = "Skew-normal PDF")

# Compare with the normal distribution with the same (xi, omega)
lines(x, dnorm(x, mean = xi, sd = omega), lwd = 2, lty = 2)

legend("topright",
       legend = c("Skew-normal", "Normal (same xi, omega)"),
       lwd = 2, lty = c(1, 2), bty = "n")

ここでは「正規分布と同じ位置・尺度をもつが、非対称性（歪み）だけが追加される」ことを視覚的に確認できます。

3.3 乱数生成

　乱数生成は rsn() を使います。生成した標本のヒストグラムに理論密度 dsn() を重ねると、分布の形が直感的に理解できます。


# ============================================================
# 2) Random sampling + histogram with theoretical PDF
# ============================================================
set.seed(1)

n <- 10000
z <- rsn(n, xi = xi, omega = omega, alpha = alpha)

hist(z, breaks = 60, freq = FALSE,
     col = "gray90", border = "gray70",
     main = "Random samples from Skew-normal",
     xlab = "x")

lines(x, dsn(x, xi = xi, omega = omega, alpha = alpha), lwd = 2)

3.4 パラメータ推定（最尤推定）

　Skew Normal の推定は、sn の selm()（Skew Elliptical Linear Model）で行うのが基本です。分布推定だけを行う場合は、z ~ 1（定数項のみ）と書きます。


# ============================================================
# 3) Parameter estimation (MLE) using selm()
# ============================================================
fit <- selm(z ~ 1, family = "SN")
summary(fit)

このスクリプトの実行結果は以下のようになります。

Call: selm(formula = z ~ 1, family = "SN")
Number of observations: 10000 
Family: SN 
Estimation method: MLE
Log-likelihood: -9551.56 
Parameter type: CP 

CP residuals:
     Min       1Q   Median       3Q      Max 
-1.85173 -0.47543 -0.08803  0.40058  2.73268 

Regression coefficients
      estimate   std.err   z-ratio Pr{>|z|}
mean 7.670e-01 6.496e-03 1.181e+02        0

Parameters of the SEC random component
       estimate std.err
s.d.     0.6529   0.005
gamma1   0.6451   0.017

この結果は、Skew Normal 分布の「位置・ばらつき・歪み」が同時に推定されたことを示しています。selm() の出力では、推定結果が CP（Centered Parameters）表現でまとめられており、これは「分布の平均・標準偏差・歪度」を直接解釈できる形になっています。
　まず、Regression coefficients に表示されている mean は、分布の平均値そのものの推定結果です。この例では平均が約 0.77 と推定されており、データが 0 を中心とする対称分布ではないことが分かります。
　次に、Parameters of the SEC random component に表示されている s.d. は、歪みを考慮した 実効的な標準偏差を表します。また、gamma1 は 分布の歪度（skewness）であり、正の値をとっていることから、分布が右に歪んでいることが定量的に確認できます。いずれも標準誤差が小さく、推定が安定していることが分かります。
　このように selm() を用いると、Skew Normal 分布の形状を特徴づける主要な量を、正規分布では捉えられない非対称性も含めて一度に推定できる点が大きな利点です。
　selm() の推定値は内部表現（最適化の都合の表現）で保持されるため、解説では 直接パラメータ（xi, omega, alpha）として取り出すのが親切です。


# Extract estimated parameters as (xi, omega, alpha)
dp <- coef(fit, param.type = "DP")
dp
# dp["xi"], dp["omega"], dp["alpha"]

推定した分布を、ヒストグラム上に重ねて確認します。


# ============================================================
# 4) Overlay fitted PDF
# ============================================================
xi_hat    <- dp["xi"]
omega_hat <- dp["omega"]
alpha_hat <- dp["alpha"]

hist(z, breaks = 60, freq = FALSE,
     col = "gray90", border = "gray70",
     main = "MLE fit: Skew-normal",
     xlab = "x")

lines(x, dsn(x, xi = xi_hat, omega = omega_hat, alpha = alpha_hat), lwd = 2)

legend("topright",
       legend = c("Data (hist)", "Fitted skew-normal"),
       lwd = c(NA, 2), pch = c(15, NA),
       col = c("gray70", "black"),
       bty = "n")

■ おまけ

　２次元正規分布を切り取ってSkew normalを作るデモが以下のRスクリプトです。

############################################################
# Skew-normal as a selection from 2D normal (explanatory figure)
#   (X,Y) ~ N(0,I), keep X only if Y <= alpha * X
# Then X | (Y <= alpha X) has density: f(x)=2*phi(x)*Phi(alpha x)
#
# Figure layout (2x2):
#  [1,1] 2D cloud + selection boundary
#  [2,1] Gate Phi(alpha x)
#  [1,2] Marginal X and Conditional X|A : KDE + theoretical curves
#  [2,2] Theoretical curves only
#
# base R only
############################################################

## -----------------------------
## 0) Parameters
## -----------------------------
n     <- 100000
alpha <- 3
m     <- 4                 # plot range for x, y
n_scatter <- 6000          # points for 2D scatter

## -----------------------------
## 1) Generate 2D standard normal
## -----------------------------
X <- rnorm(n)
Y <- rnorm(n)

A <- (Y <= alpha * X)       # selection event
X_acc <- X[A]
Y_acc <- Y[A]

cat("Acceptance rate ~", mean(A), "\n")  # ~ 0.5

## -----------------------------
## 2) Theoretical functions
## -----------------------------
phi   <- function(x) dnorm(x)
Phi   <- function(x) pnorm(x)
f_skew <- function(x, alpha) 2 * dnorm(x) * pnorm(alpha * x)
delta <- alpha / sqrt(1 + alpha^2)
mean_theory <- delta * sqrt(2 / pi)

grid_x <- seq(-m, m+mean_theory, length.out = 1200)

## -----------------------------
## 3) Density estimates (KDE)
## -----------------------------
# Use a common range for fair comparison
dens_X     <- density(X,     from = -m, to = m, n = 1024)
dens_X_acc <- density(X_acc, from = -m, to = m, n = 1024)

## -----------------------------
## 4) Plot settings
## -----------------------------
op <- par(no.readonly = TRUE)
on.exit(par(op), add = TRUE)

par(mfcol = c(2, 2), mar = c(4, 4, 3, 1),xaxs="i",pty="s")

## =========================================================
## [1,1] 2D cloud + selection half-plane
## =========================================================
idx <- sample.int(n, min(n, n_scatter))
idx_acc <- idx[A[idx]]

plot(X[idx], Y[idx],
     pch = 16, cex = 0.4, col = "gray75",
     xlim = c(-m, m), ylim = c(-m, m),
     xlab = "X", ylab = "Y",
     main = "2D Normal + Selection")

points(X[idx_acc], Y[idx_acc], pch = 16, cex = 0.5, col = "pink")
# boundary y = alpha x
abline(a = 0, b = alpha, lwd = 2,col=4)

legend("topleft",
       legend = c("All points", "Accepted", expression(y==alpha*x)),
       pch = c(16, 16, NA), lty = c(NA, NA, 1),lwd=c(NA,NA,2),
       col = c("gray75", "pink", 4),
       bty = "n", cex = 0.85)

## =========================================================
## [2,1] Gate Phi(alpha x)
## =========================================================
plot(grid_x, Phi(alpha * grid_x),
     type = "l", lwd = 2,
     xlim = c(-m, m), ylim = c(0, 1),
     xlab = "x", ylab = expression(Phi(alpha*x)),col=4,
     main = expression(paste("Gate  ", Phi(alpha*x))))
abline(h = 0.5, lty = 2, col = "gray")
abline(v = 0,   lty = 2, col = "gray")

## =========================================================
## [1,2] Marginal vs Conditional (KDE + theory)
## =========================================================
yl <- range(
  dens_X$y, dens_X_acc$y,
  phi(grid_x), f_skew(grid_x, alpha)
)*1.3

plot(dens_X$x, dens_X$y,
     type = "l", lwd = 2, col = NA,
     xlim = c(-m, m), ylim = yl,
     xlab = "x", ylab = "Density",
     main = "Monte Carlo vs Theoretical Distributions")

polygon(c(dens_X$x,dens_X$x[1]), c(dens_X$y,0),col="gray75",border=NA)

# KDE (accepted)
polygon(c(dens_X_acc$x,dens_X_acc$x[1]), c(dens_X_acc$y,0),col=rgb(1,0.4,0.4,alpha=0.3),border=NA)

# Theory curves
lines(grid_x, phi(grid_x), lwd = 2, lty = 2, col = "gray40")
lines(grid_x, f_skew(grid_x, alpha), lwd = 2, lty = 2, col = 2)

legend("topleft",
       legend = c("Density estimetion of X",
                  "Density estimetion of X | (Y <= alpha X)",
                  "Theory: phi(x)",
                  expression("Theory: 2 phi(x) Phi(alpha x)")),
       col  = c("gray75", rgb(1,0.4,0.4,alpha=0.3), "gray40", 2),
       lwd  = c(NA, NA, 2, 2),pch=c(15,15,NA,NA),
       lty  = c(NA, NA, 2, 2),
       bty  = "n", cex = 0.85)

## =========================================================
## [2,2] Theory curves only
## =========================================================


plot(grid_x-mean_theory, f_skew(grid_x, alpha),
     type = "l", lwd = 2, col = 2,
     xlim = c(-m, m),
     ylim = range(phi(grid_x), f_skew(grid_x, alpha))*1.3,
     xlab = "x", ylab = "Density",
     main = "Zero-mean comparison")

lines(grid_x, phi(grid_x), lwd = 2, lty = 2, col = "gray40")

legend("topleft",
       legend = c(expression("Skew-normal with zero mean"),
                  expression("Normal")),
       col = c(2, "gray40"),
       lwd = 2, lty = c(1, 2),
       bty = "n", cex = 0.85)

par(op)

4. まとめ

　今回は、非対称分布の例としてskew normalを紹介しました。紹介した理由は、正規性検定や正規性の可視化をためすため、非対称分布を導入する論文を書くときに先行研究としてskew normalを紹介するため、自分のデータにskew normalを適用してみるためです。　Skew normalのwikipediaの解説を読んだり、式を見つめたりするだけでは、理解は深まりませんので、Ｒで自分で体験してみてください。

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