📂確率分布論

パレート分布

定義 ¹

スケールパラメータ$x_{0} > 0$とシェイプパラメータ$\alpha > 0$の場合、以下の確率関数をパレート分布、パワー法則、またはスケールフリー分布という。

1. 連続型：定数$\displaystyle \int_{x_{0}}^{\infty} p(x) dx = 1$を満たすための定数$C$について $$p(x) = C x^{-\alpha} \qquad , x > x_{0}$$
2. 離散型：リーマンのゼータ関数$\zeta$について $$p_{k} = {{ 1 } \over { \zeta (\alpha) }} k^{-\alpha} \qquad , k \in \mathbb{N}$$

基本的な性質

• [1] モーメント生成関数：パレート分布のモーメント生成関数は存在しない。
• [2] 平均と分散：$X \sim \text{Pareto} \left( x_{0}, \alpha \right)$の場合 $$\begin{align*} E (X) =& {{ \alpha - 1 } \over { \alpha - 2 }} x_{0} & , \alpha > 2 \\ \Var (X) =& {{ (\alpha - 1) } \over { \left( \alpha -2 \right)^{2} (\alpha - 3) }} x_{0}^{2} & , \alpha > 3 \end{align*}$$

定理

• [a] スケールフリー性：パレート分布は唯一のスケールフリー分布だ。つまり、すべての$b$についてある定数$\alpha$が存在し、次が成立する。 $$p(bx) = g(b) p(x) \implies p(x) = p(1) x^{-\alpha}$$
• [b] $k$次のモーメント：$0 < k < \alpha - 1$の場合、$k$次のモーメントが存在し $$E X^{k} = {{ \alpha - 1 } \over { \alpha - 1 - k }} x_{0}^{k}$$

説明

パレート分布はこの実際の世界に蔓延っている不平等を説明する代表的な分布で、次の概念と非常に密接な関係がある。

• ヒープスの法則とジップの法則：単語の頻度に関する経験則。
• 本の販売量
• 通信量
• 地震の強さ
• クレーターの直径
• 財産
• 引用数

確率密度関数の形を見ると、シェイプ$\alpha$が大きいほど不平等が激しくなることが直感的に理解できる。経済的に言えば、富裕層はお金が無限にあり、貧しい人がたくさんいるということだ。

パレート分布がスケールフリー性を持つという話は、文字通りスケールがないということだ。例えば、ポアソン分布に従う2つの確率変数のパラメータが$\lambda_{1} = 10$、$\lambda_{2} = 1000$であれば、見る場所によって大きな違いがあるが、パレート分布はどこを見ても本質的に違いがないということだ。数学的には$b$がどのような値で与えられても結論が同じであることに相当する。

証明

[1]

確率変数のモーメント生成関数が存在することは、すべての$k \in \mathbb{N}$に対して$k$次のモーメントが存在することを意味する。しかし、定理[2]で示されたように、パレート分布の$1$次のモーメントは$\alpha > 1$の場合にのみ存在するので、モーメント生成関数は存在しえない。

■

[2]

戦略：モーメントの公式[b]を利用する。

---

$$\begin{align*} EX^{1} =& {{ \alpha - 1 } \over { \alpha - 1 - 1 }} x_{0}^{1} \\ =& {{ \alpha - 1 } \over { \alpha - 2 }} x_{0}^{1} \end{align*}$$ であり、$\displaystyle EX^{2} = {{ \alpha - 1 } \over { \alpha - 3 }} x_{0}^{2}$なので $$\begin{align*} \Var X =& {{ \alpha - 1 } \over { \alpha - 3 }} x_{0}^{2} - \left[ {{ \alpha - 1 } \over { \alpha - 2 }} x_{0}^{1} \right]^{2} \\ =& \left[ {{ 1 } \over { \alpha - 3 }} - {{ \alpha - 1 } \over { \left( \alpha - 2 \right)^{2} }} \right] (\alpha - 1) x_{0}^{2} \\ =& \left[ \alpha^{2} - 4 \alpha + 4 - \alpha^{2} + 4 \alpha - 3 \right] {{ (\alpha - 1) } \over { (\alpha - 3) \left( \alpha -2 \right)^{2} }} x_{0}^{2} \\ =& {{ (\alpha - 1) } \over { \left( \alpha -2 \right)^{2} (\alpha - 3) }} x_{0}^{2} \end{align*}$$

■

[a]

すべての$b$に対してある関数$g$が存在し $$p(bx) = g(b) p(x)$$ が成立すると仮定する。$x = 1$を代入してみると$p(b) = g(b) p(1)$であるため、$g(b) = p(b) / p(1)$であり $$p(bx) = {{ p(b) p(x) } \over { p(1) }}$$ $b$に対して微分してみると $$x p '(bx) = {{ p ' (b) p(x) } \over { p(1) }}$$ $b=1$を代入してみると対数関数の微分法を用いたトリックにより² $$\begin{align*} & x p '(x) = {{ p ' (1) p(x) } \over { p(1) }} \\ \implies & {{ p '(x) } \over { p(x) }} = {{ p '(1) } \over { p(1) }} \cdot {{ 1 } \over { x }} \\ \implies & {{ d \log p(x) } \over { dx }} = {{ p '(1) } \over { p(1) }} \cdot {{ 1 } \over { x }} \\ \implies & d \log p(x) = {{ p '(1) } \over { p(1) }} {{ 1 } \over { x }} dx \end{align*}$$ これは、ある定数$\text{constant}$に対する単純な分離可能な1階微分方程式で、次を得る。 $$\log p(x) = {{ p '(1) } \over { p(1) }} \log x + \text{constant}$$ $x = 1$を代入してみると$\text{constant} = \log p(1)$であることがわかる。$\displaystyle \alpha := - {{ p '(1) } \over { p(1) }}$と定義すると、求めていた次の式を得る。 $$\begin{align*} & \log p(x) = - \alpha \log x + \log p(1) \\ \implies & \log p(x) = \log x^{-\alpha} + \log p(1) \\ \implies & \log p(x) = \log x^{-\alpha} p(1) \\ \implies & p(x) = p(1) x^{-\alpha} \end{align*}$$

■

[b]

$0 < \alpha -1$であるので、$\displaystyle \int_{x_{0}}^{\infty} C x^{-\alpha} dx = 1$から$C = \left( \alpha - 1 \right) x_{0}^{\alpha - 1}$を得る。したがって、 $$\begin{align*} E X^{k} =& \int_{x_{0}}^{\infty} x^{k} C x^{-\alpha} dx \\ =& C \int_{x_{0}}^{\infty} x^{k-\alpha} dx \\ =& \left( \alpha - 1 \right) x_{0}^{\alpha - 1} \left[ {{ 1 } \over { k - \alpha + 1 }} x^{k - \alpha + 1} \right]_{x_{0}}^{\infty} \\ =& \left( \alpha - 1 \right) x_{0}^{\alpha - 1} \left( 0 - {{ 1 } \over { k - \alpha + 1 }} x_{0}^{k - \alpha + 1} \right) \\ =& {{ \alpha - 1 } \over { \alpha - 1 - k }} x_{0}^{k} \end{align*}$$

■

可視化

次に、パレート分布の確率密度関数をGIFで表示するJuliaのコードです。

@time using LaTeXStrings
@time using Distributions
@time using Plots

cd(@__DIR__)

x = 1:0.1:10
A = collect(0.5:0.01:3.5); append!(A, reverse(A))

animation = @animate for α ∈ A
    plot(x, pdf.(Pareto(α), x),
     color = :black,
     label = "α = $(round(α, digits = 2))", size = (400,300))
    xlims!(0,5); ylims!(0,4); title!(L"\mathrm{pdf\,of\,Pareto}(\alpha)")
end
gif(animation, "pdf.gif")