「R」に関する備忘録

in [数ナビの部屋]

「R」に関する諸操作を備忘録的に取りまとめました.

統計処理ソフト「R」を活用して統計解析を究めよう!

  • モバイル利用(Android)でのメニュー選択は、 SiteMapを利用するか、 または「長押し」から「新しいタブを開く」を選択してください。
■ 「R」に関する備忘録
[御案内] フリーの本格的な統計処理ソフトに「R」があります. ここでは,「R」のパッケージを取り扱うに至るまでの 「R」の利用法について備忘録的に取りまとめました. 「R Studio」は使用していません.



パッケージ

パッケージのインストール
「R」では,標準で備わっている機能の他に, パッケージを読み込ませることでいろいろな機能を利用することができる. 標準で登録されているパッケージは「library()」とすると表示される.

ただし,表示されるすべてのパッケージをそのまま利用できるわけではない. その時点で利用できるパッケージは「search()」により表示される. 下図は,3列で表示される結果の2列目までを切り取ったものである.

起動直後では「stats」「graphics」「grDevices」「utils」 「datasets」「methods」「base」が読み込まれている. 特殊な使い方をしなければ, これだけで普通の統計解析やグラフ表示はカバーされるのでは ないかと思われる.「library()」で表示されても「search()」で 表示されないパッケージを読み込ませるには, 画面の上部メニューの「パッケージ」から「パッケージの読み込み」を 選択する.読み込み可能なパッケージ名が表示されるので, それを選択する.読み込ませたいパッケージ名が分かっているときは, 直接「library(パッケージ名)」を実行してもよい.

上記の「パッケージの読み込み」で表示されなくても, 「R」で利用できるパッケージとしては膨大な量のパッケージがある. どのようなパッケージを利用できるかは, 上部メニューの「パッケージ」から「パッケージのインストール」を選択する. 最初に,ダウンロードできるサーバーのリストが表示されるので, 希望するサーバーを選択する. 普通に「Japan(Tokyo)[https]」を選択すればよい. そのサーバーに接続するには多少の時間がかかるので, 10秒前後はじっと待つ.接続すると, ダウンロード可能なパッケージの膨大なリストが表示されるので, そこからパッケージを選択する. なお,当然ながら, 自分が必要とする機能を持つパッケージはどのようなパッケージであるかは, 事前に調べていることが前提となる.

「R」で標準で利用できる確率分布に「ベキ分布(パレート分布)」は 登録されていないので,ここでは「Pareto」を選択する. なお,「ベキ分布」のMaximaでの解説は 「こちら」を参照されたい.

選択して「OK」を押すと,パッケージの保存場所が表示される. 隠しフォルダーである「AppData」に保存されるようである.

そのパッケージは自動的に展開されてチェックサムも行われ, 「R」の作業フォルダーに保存される. その作業フォルダーには「win-library」というフォルダーがあり, その下部フォルダーにバージョン別のフォルダーがある. 使用している「R」のバージョンは「3.6.2」なので「3.6」の フォルダーに進むと読み込まれたパッケージ別のフォルダーがある. そこで,「Pareto」のフォルダーに進むと「doc」フォルダーがあり, そこに「Pareto」パッケージの解説が収められている. つまり,パッケージにもよると思われるが,「Pareto」パッケージの場合の 解説は, [Rの作業フォルダー]>[win-library]>[3.6]>[Pareto]>[doc] の箇所にある「Pareto.html」である. 下図は,その冒頭部分である.日付が「2020-02-13」となっているので, 最近更新されたドキュメントであることが分かる (この文面は「2020-03-22」に書いている).

ただし,パッケージをインストールしても, それは読み込まれるだけで利用可能になるわけではない. 利用するには,改めて「library(Pareto)」とする.

ヘルプファイルがあれば「help(Pareto)」で表示される. この場合は無いようなので,メッセージにしたがい「??Pareto」とすると, 関数名に「Pareto」を含む「R」のパッケージのリストが表示される.

これを見ると,「Pareto」という名前を含む関数を持つパッケージは, 「Pareto」の他に「actuar」「lmon」「ParetoPosStable」「rmutil」「VGAM」が あることが分かる.右側の解説を読むと,パレート分布に関する 「密度関数(dPareto)」「分布関数(pPareto)」「分位関数(qPareto)」 「乱数(rPareto)」をすべて含んでいるのは「Pareto」パッケージだけである ことが分かる.左側をクリックすると, その関数の引数や具体的な使用例が表示される. 下図は「Pareto:dPareto」をクリックした場合である.


▲戻る(トップメニュー)

ベキ分布(パレート分布) [script]
以下は,「Pareto.html」をもとに, このパッケージに登録されている関数のうち,下記の関数の使い方について解説する.

■確率密度関数(dPareto)
ベキ分布(パレート分布) Pareto(t,α) の累積分布関数 \(\small F_{t,\alpha}(x)\) は, 次の式で表される. \[\small F_{t,\alpha}=\begin{cases} \quad 0 & (x\leq t)\\ \displaystyle 1-\left(\frac{t}{x}\right)^{\alpha}& (x>t) \end{cases} \] これを微分すれば確率密度関数が得られる. \[\small \frac{d}{dx}F_{t,\alpha}(x) =\begin{cases} \quad 0 & (x\leq t)\\ \displaystyle \frac{\alpha t^{\alpha}}{x^{\alpha+1}} & (x>t) \end{cases}\] 「パレート分布」といえば,通常は上記の関数で定義される確率分布を指す. さらに専門的になると,この分布は「パレートI型」,「ヨーロピアンパレート」, あるいは「1パラメータパレート」とも呼ばれるようであるが, 以後,本サイトでは「ベキ分布」と呼ぶ. また,「Pareto」パッケージをインストールして 「library(Pareto)」を実行済みとする.

「R」での関数名は, 確率密度関数は「dPareto」,累積分布関数は「pPareto」である. Maximaでの確率密度関数は \(\small ab^{a}/x^{a+1}~(x\geq b)\) で考えて主に \(\small b=1\) の場合を扱ったので, ここでも \(\small t=1\) の場合を考える.

「dPareto」の書式は 「dPareto(x, t, \(\small \alpha\), truncation)」である. 「truncation」は省略可能で, デフォルトでは「truncation=NULL」になっている. 「truncation=s」とすると,\(\small t\lt x\le s\) の部分のグラフが描画され, \(\small x\gt s\) のときは「0」となる. また,「x」に点列を指定するベクトルを与えると, 「plot」で描画することができる. なお,開始点がMaximaでは \(\small x\geq b\) であるが, 「R」では \(\small x>t\) であるので注意が必要である. 範囲を \(\small x\geq t\) で指定すると,値が「0」からのグラフとなる.

plot(dPareto(1:20,1,1),type="o")
\(\small 1,2,3,\cdots,20\) のときの「dPareto」の値が 点で配置され,「type="o"」を指定したので線で繋いで描画される. \(\small x\le t\) では「dPareto=0」なので,グラフは横軸から始まる.

plot(dPareto(1.01:20,1,1),type="o")
点列を 1.01:20 としたので,実際には \(\small 1.01, 2.01, \ldots,20.01\) のときの値で描画される.

curve(dPareto(x,1,1),xlim=c(1,20))
範囲を \(\small [1,20]\) としたので,値が「0」からのグラフとなる.

curve(dPareto(x,1,1,5),xlim=c(1.01,20))
「truncation=5」を指定したので,\(\small x\gt 5\) の部分は「0」となる.

以下では,指数 \(\small \alpha\) の違いをみるため, \(\small \alpha=0.5(赤), 1(黒), 1.5(緑)\) の場合のグラフである.

以下は,同じグラフで縦軸と横軸を対数軸としたものである. ベキ分布は,両対数グラフでは直線になるのが特徴的である.

▲戻る(トップメニュー)

■累積分布関数(pPareto)
ベキ分布の累積分布関数 \[\small F_{t,\alpha}(x) =\begin{cases} \quad 0 & (x\leq t)\\ \displaystyle 1-\left(\frac{t}{x}\right)^{\alpha}& (x>t) \end{cases}\] の「R」での名称は「pPareto」であり,書式は 「pPareto(x, t, \(\small \alpha\), truncation)」である. 使い方は「dPareto」と同様である. 累積分布関数の初期値は「0」であるので, 端点を気にする必要はない.

plot(pPareto(1:20,1,1),type="o")

以下は,指数 \(\small \alpha\) の違いをみるため, \(\small \alpha=0.5(赤), 1(黒), 2(緑)\) の場合のグラフである. 横軸を対数軸とした.

▲戻る(トップメニュー)

■分位関数(qPareto)
分位関数「qPareto」は,累積分布関数「pPareto」の逆関数である. つまり,「\(\small {\rm qPareto}\,(q,t,\alpha)=x\)」とすると, これは「\(\small {\rm pPareto}\,(x,t,\alpha)=q\)」,すなわち \(\small X\) を Pareto(\(\small t,\alpha\)) にしたがう確率変数とするとき, \(\small {\rm P}(X\le x)=q\) ということである.

上図で「q <- 1:5/5」は,ベクトル c(1,2,3,4,5) を「5」で割って, c(0.2, 0.4, 0.6, 0.8, 1) というベクトルを「q」とした. 各点での分位数が得られる.念のため,c(2.5, 5.0) について pParetoの値を確かめた.

下図は,qPareto(x,1,1)のグラフである. 縦軸を対数軸とした.

curve(qPareto(x,1,1),xlim=c(0,1),log="y")

▲戻る(トップメニュー)

■乱数(rPareto)
「rPareto」を利用すると, ベキ分布 Pareto(\(\small t,\alpha\)) にしたがう 乱数を発生させることができる.

下図では,Pareto(1,1) にしたがう乱数1000個からなるベクトルを 「pare2」とした.「summary」を調べて左側だけを切り出した. 第3四分位数(3rd Qu.)は 3.984,最大値は 1062.54である. 次に,区間 [1,1063] を1刻みにしてヒストグラムを作成し, 区間 [1.01,15] の範囲で確率密度関数「dPareto(x,1,1)」の グラフを「add=T」として重ね合わせた.


▲戻る(トップメニュー)

■標本平均の分布
dPareto(x,1,1) は \(\small f(x)=\frac1{x^2}\) ということである. 簡単な計算 から分かるように,この確率密度関数は平均も分散も存在しない. 中心極限定理によれば,平均や分散が存在するような確率分布では, 標本数を増やしていくと標本平均の分布は正規分布で近似される.

ここでは,平均や分散が存在しない「ベキ分布」にしたがう 乱数を「rPareto(n,1,1)」により発生させて, その標本平均がどのような分布になるかを調べる. 最初に,n個の標本平均からなる要素数N個のベクトルを生成する関数を定義する.

次に,この関数を利用して,n=2,5,10,100,1000,10000 の場合の 標本平均を要素とするベクトルを生成する.

1万個の平均を千回計算するのはかなりの時間がかかりそうに見えるが, ほぼ一瞬で生成された.ちなみに, Maximaでこの計算をさせると,かなりの時間を要する. さすがは統計処理ソフト「R」と実感させられた. 次に,これらを列ベクトルにしてまとめて,その「summary」を 表示させる.

次に,これらのデータのヒストグラムを折れ線で同時に描画するため, それぞれの最大値の値に応じて1刻みとしたヒストグラムのデータを取得する. 念のため,hp2 の列名を表示させた.

以上のデータをもとに「ヒストグラム(2)」 で述べた折れ線で描画する.色は, n=2(黒), 5(赤),10(緑), 100(青),1000(水色),10000(紫) とした.

グラフを区別するために低水準作図関数の「text」を利用すると, 「text(x,y,文字)」の形で文字を書き入れることができる. ただし,座標の指定には多少の試行錯誤が必要である. こみ入ったグラフをそのつど書き直すのは面倒なので, スクリプトに書いておいて「ctrl+R」で一気に描画するのがよい.

これを見ると,標本数が増えるとグラフは右側にずれていくが, グラフは同じような形に近づいているように思われる. 念のため n=100,000 の場合も計算すると, さすがに数秒の時間を要したが同じようなグラフが n=10,000 のグラフの 右側に表示された(図は略).

実は,ベキ分布にしたがう確率変数の平均は「安定分布」に 近づいていくことが証明されている(参照). また,Pareto(1,2/3)とPareto(1,3/2)の場合は, 標本平均の収束先の確率密度関数を特殊関数を用いて 表すことができる.それを実際に行ったのが, 次のパッケージ「gsl」での解説である.


▲戻る(トップメニュー)

GNU Science Library(gsl)  [script]

「R」のパッケージに「gsl」がある. これは「GNU Science Library(GSL)」と呼ばれるもので, 他方面にわたる数学関数を取りまとめたC言語のライブラリーである. 「gsl」は,そのライブラリーに含まれる幾つかの関数を 「R」でも利用できるようにしたものである。

■「GSL」とは
「GSL」ではどのような関数を利用できるかは下記を参照するとよい.
  1. GNU科学技術計算ライブラリー,リファレンス・マニュアル[578頁]
  2. Package 'gsl'

[1]は「GSL」の解説書(和文)であり膨大な関数群が登録されている. [2]は「R」での利用例が紹介されている(英文). なお,「GSL」に登録されている関数は, C言語に堪能であれば「gnuplot」で 利用することもできる.

■合流型超幾何関数
ここでは,「ベキ分布」での解析で必要となることから 合流型超幾何関数の使用例を述べる. 合流型超幾何関数は合流型超幾何微分方程式 \[\small zy''+(b-z)y'-ay=0\] の解として定義される. この微分方程式の一般解は \[\small y=c_1M(a,b,z)+c_2U(a,b,z)\] で表され, \(\small M(a,b,z)\) は第1種合流型超幾何関数, \(\small U(a,b,z)\) は第2種合流型超幾何関数 と呼ばれ,次式で定義される. (参照1参照2). \[\small \begin{align*} M(a,b,z) &=1+\frac{a}{b}z+\frac{a(a+1)}{b(b+1)}\frac{z^2}{2!}\\ &+\cdots=\sum_{n=0}^{\infty}\frac{(a)_n}{(b)_n}\frac{z^n}{n!}\\ U(a,b,z) &=\frac{\Gamma(1-b)}{\Gamma(a-b+1)}M(a,b,z)\\ &+\frac{\Gamma(b-1)}{\Gamma(a)}z^{1-b}\\ &\quad \times M(a-b+1,2-b,z) \end{align*}\] なお,\(\small (a)_n=a(a+1)(a+2)\cdots (a+n-1)\) であり, \(\small \Gamma(z)\) はガンマ関数である.

この第2種合流型超幾何関数は,「gsl」では「hyperg_U」として登録されている. さらに,この関数を用いると,「(第2種)ホイッタカー(Whittaker)関数」 なるものが次式で定義される.この関数は,次の項目で述べる ベキ分布にしたがう確率変数の和の分布を考えるときに必要となる.

\[\small \begin{align*} W_{k,m}(z)&=e^{-\frac{z}{2}}z^{m+\frac12}{}\\ &\times U(m-k+\frac12,1+2m,z)\end{align*}\] たとえば,\(\small w(2,1/2,x)\) のグラフは下図となる. WolframAlphaでも同じグラフが描画されている.


ここで,「abline」は低水準作図関数であり, 「abline(0,0)」で x 軸を,「abline(v=0)」で y 軸が描かれている. 「abline(a,b)」で直線 \(\small y=a+bx\) が描かれ, 「abline(v=x)」で指定した x 座標の箇所に垂線が引かれる. 「abline(h=y)」とすると,指定した y 座標の箇所に水平線が引かれる. したがって,「abline(0,0)」のかわりに「abline(h=0)」としてもよい.
▲戻る(トップメニュー)

■ベキ分布の和
このホイッタカー関数を利用すると, ベキ分布 \(\small f(x)=\frac{\alpha}{x^{\alpha+1}}~(x\gt 1, \alpha\gt 0)\) にしたがう確率変数の和の分布について考えることができる.

ベキ分布の確率密度関数は単純であるが, この分布にしたがう確率変数の和の分布は簡単な式で表すことはできない. この分布は指数 \(\small \alpha\) の値によっては 平均や分散が存在するとは限らないため,中心極限定理は適用できない. しかし,分散が存在しない確率分布の場合の中心極限定理に対応するものとして, 「一般化中心極限定理」がある.それによると, ベキ分布にしたがう互いに独立な確率変数の和を \(\small S_n\) とすると, 定数 \(\small a_n, b_n\) に対して \[\small \lim_{n\to\infty} {\rm P} \left(\frac{S_n-b_n}{a_n}\leq x\right) =F_{\alpha}(x)\] を満たすような安定分布 (Rでの解説)が存在することが証明されている (参照, p.12(19)) . 右辺は,指数 \(\small \alpha\) を持つ安定分布(Maximaでの解説)の累積分布関数である.(実際には,他に 3つのパラメータを持つ.) ここで,定数 \(\small a_n, b_n\) は \(\small a_n=n^{1/{\alpha}}C_{\alpha}\) であり,\(\small C_{\alpha}\) は次式で表される. \[\small C_{\alpha} =\begin{cases} \small \left(\Gamma(1-\alpha)\cos\left(\frac{\pi \alpha}{2}\right)\right)^{1/{\alpha}} & (\alpha\ne 1) \\ \displaystyle \quad \frac{\pi}{2} &(\alpha=1) \end{cases}\] \(\small b_n\) は,\(\small 0\lt \alpha\lt 1,~ \alpha=1,~ 1\lt\alpha \lt 2\) に応じて,それぞれ \[\small b_n=\begin{cases}\small\quad 0\\ \displaystyle\small n\log{n}+n\left(1-C-\log\frac2{\pi}\right)\\ \displaystyle\small \frac{n\alpha}{\alpha-1} \end{cases}\] ここで,\(\small b_n\) の2行目は近似式であり,\(\small C\) はオイラー定数 (\(\small 0.5772156649\cdots\))である.

通常の統計の授業で安定分布に触れられることはないが, 現実社会のいろいろな現象を解析する上で安定分布は極めて重要な確率分布である. しかし,その確率密度関数を初等関数で書き表すことができるのは, 正規分布,コーシー分布,そしてレビ分布に限られ, それぞれ \(\small \alpha=2, 1, 0.5\) の安定分布に対応する (後半の2つの分布は平均も分散も存在しない). しかし,特殊関数であるホイッタカー関数を利用すると, \(\small \alpha=\frac23, \frac32\) の場合も式で書き表すことができる. それぞれの確率密度関数を \(\small f_{2/3}(x), f_{3/2}(x)\) とすると, 次の式で表される(参照,p.10(15〜18)).

\(\small f_{2/3}(x)~(x\gt 0)\) は
\[\small \sqrt{\frac{3}{\pi}}\frac1{x}\exp\left(\small -\frac{16}{27x^2}\right) W_{1/2,1/6}\left(\small \frac{32}{27x^2}\right)\] \(\small f_{3/2}(x)\) は,\(\small x\lt 0, x\gt 0\) の場合に, それぞれ \[ \begin{cases}\small -\sqrt{\frac{3}{\pi}}\frac1{x}\exp\left(\frac{x^3}{27}\right) W_{1/2,1/6}\left(-\frac{2x^3}{27}\right)\\ \small \frac1{2\sqrt{3\pi}}\frac1{x}\exp\left(\frac{x^3}{27}\right) W_{-1/2,1/6}\left(\frac{2x^3}{27}\right) \end{cases} \] 以下では,\(\small \alpha=\frac23, \frac32\) の場合の ベキ分布 \(\small f(x)=\frac{\alpha}{x^{\alpha+1}}\) にしたがう乱数を発生させて,このことを実際に確かめてみる.

[1] \(\small \alpha=\frac23\) の場合
最初に,確率密度関数を定義して,同じ指数をもつベキ分布と比較する. \(\small x\gt 2\) ではほぼ重なっているのが分かる.


次に,定数 \(\small a_n, b_n\) を定義する.

乱数 rPareto(1,2/3) を発生させて \(\small n\) 個の平均からなる 要素数 1000 のベクトルを生成する関数を定義する.

その関数を利用して,\(\small n=2, 5, 10, 100\) の平均からなる ベクトルを生成し,それを列ベクトルにしてまとめて「summary」を出力させる.

生成したベクトルから,1刻みにした場合のヒストグラム情報を取得する. 最小値と最大値の中に全データが含まれるように指定しないとエラーが表示される.

それをもとに折れ線グラフを描画すると,次のようになる. 色は,2個(黒),5個(赤),10個(緑),そして100個(青)である. 平均を取る個数が多くなるにつれてなだらかになっていく.

今度は,これを標準化して \(\small (S_n-b_n)/a_n\) について考える. 平均を取っているので,\(\small (\bar{X}_n-b_n/n)/(a_n/n)\) を考えることになる. すでに定義済みの係数を用いて標準化する関数を定義して, 生成した平均からなるベクトルを標準化する.

次に,列ベクトルにまとめて「summary」を実行して最小値と最大値を確認する.

以上の値をもとに,0.5刻みのヒストグラム情報を取得して,折れ線グラフを描画する. 色分けは,2個(黒),5個(赤),10個(青),100個(緑),そして 確率密度曲線 \(\small f_{2/3}(x)\) (シアン)である.個数が増えるにつれ, 密度曲線に近づいていくのが分かる.

下図は,1000個に増やして棒グラフで描画したものである. この場合の最大値は 6785.806 である.

▲戻る(トップメニュー)

[2] \(\small \alpha=\frac32\) の場合
確率密度関数を定義して,同じ指数をもつベキ分布と比較する. 正負で異なる関数で定義されるので,0の近くの双方の値を調べることで 「0.1975162」を決定した.同じ指数のベキ分布 \(\small f(x)=\frac{3/2}{x^{3/2+1}}\) のグラフも描画した.[1]の場合と比較すると, ある程度の重なりを見るには値がかなり大きい必要がある.


下図は,\(\small n=10, 100, 1000\) の場合の, 標本平均の 0.5刻みとしたヒストグラム(折れ線)である. この指数の場合は平均が存在するので, 個数を増やすほど平均の近くに寄ってくることになる. 色分けは,10個(黒),100個(赤),そして1000個(緑)である.

下図は,同じ個数の場合を [1] と同様の係数で標準化したときのグラフである. [1] と比べると収束が遅いことが分かる.

下図は,\(\small n=10,000\) の場合について標準化して棒グラフで 表示したものである.さらに密度関数に近づいていることが分かる.

★以上により,指数 \(\small \alpha=\frac23, \frac32\) をもつベキ分布の 標本平均が,標本数を増やしていくと同じ指数を持つ安定分布に近づいていくことを 確認することができた.
▲戻る(トップメニュー)

安定分布  [script]


■安定分布とは
通常の統計の教科書では平均や分散が存在するような確率分布しか扱われないが, 実際の現象では「ベキ分布」のように 平均や分散が存在しない確率分布がある. そのような確率分布では中心極限定理は適用できないので, いろいろな統計量を正規分布に帰着させることができない. 平均や分散が存在しない場合には「一般化中心極限定理」があり, 「安定分布」と呼ばれる確率分布が中心になる.

安定分布は,確率変数 \(\small X\) の独立なコピーを \(\small X_1, X_2,\ldots,X_n\) とするとき, \[\small X_1+X_2+\cdots+X_n\overset{d}{=}a_nX+b_n\] を満たすような定数 \(\small a_n, b_n\) が 存在するような確率分布として定義される. 記号「\(\small \overset{d}{=}\)」は,両辺の分布が一致することを示す.

正規分布は,この性質を満たすので安定分布である.実際, 確率変数 \(\small X\) が正規分布 N\(\small (\mu,\sigma^2)\) にしたがえば, 標準化した \(\small Z=(X-\mu)/\sigma\) は \(\small X=\sigma Z+\mu\) と表され, \(\small Z\) は標準正規分布 N(0,1)にしたがう. さらに,大きさ \(\small n\) の標本平均 \(\small \bar{X}\) も 正規分布 N(\(\small \mu,\sigma^2/n\)) にしたがうので, それを標準化した \(\small (\bar{X}-\mu)/(\sigma/\sqrt{n})\) も 標準正規分布 N(0,1) にしたがう.分布が一致するので \[\small \frac{\bar{X}-\mu}{\sigma/\sqrt{n}} \overset{d}{=}Z\] と表せる. \(\small \sigma Z=X-\mu\) を代入して 整理すると \[\begin{align*}\small \small X_1+X_2+&\small \cdots+X_n\\ &\small \overset{d}{=}\sqrt{n}X+(n-\sqrt{n})\mu \end{align*}\] と表される. これは,同一の確率分布にしたがう独立な確率変数の和が 同じ分布にしたがう確率変数の1次式の分布と一致することを示している. 安定分布は,この関係を平均や分散が存在しないような確率分布に 拡張して考えて得られたものである.

安定分布は聞き慣れない確率分布であるが,平均や分散が存在する正規分布や, それらが存在しないコーシー分布をも含む広い範疇の確率分布である. この確率分布は,確率密度関数をフーリエ変換した式で特徴づけられる. それをフーリエ逆変換して確率密度関数を初等関数で表すことができるのは, 正規分布,コーシー分布,そしてレビ分布に限られる. 詳細は「こちら」(Maxima利用)を参照されたい.


▲戻る(トップメニュー)

■安定分布のパッケージ「fBasics」
安定分布を扱うパッケージとして,「R」には「fBasics」 が登録されている(ドキュメント[pdf]). しかし,このパッケージを読み込むだけでは不十分で, 他に「stabledist」「mgcv」「tcltk」も読み込む必要がある. さらに,ヒストグラムを描画しやすい「truehist」を使うために「MASS」も 読み込んでおいた方が良い. 以下では,それらを読み込み済みとする.

「fBasics」でWeb検索すると, 富士通のパソコンに組み込まれていた「FBASIC」に関するものが 多数検索されるが, 「fBasics」は経済関係(Finance)のデータを分析するための基本パッケージであり, 膨大なデータを解析するための多数の関数が登録されており, 列ごとの分散や標準偏差などをまとめて求めることもできる(colVars,colSds). さらに,一連の関数群の中に安定分布に関する下記のような関数がある.

  • dstable:確率密度関数
  • pstable:累積分布関数
  • qstable:分位関数
  • rstable:乱数生成
上記は,実際には「stabledist」に登録されている関数である. これらを使用するだけであれば,「stabledist」を読み込むだけでよい. 安定分布は4つのパラメータを持ち「dstable」であれば, 「dstable(x,\(\small \alpha, \beta, \gamma, \delta\))」の書式である. それぞれのパラメータの意味は,次のようなものになる.
  • 特性指数「\(\small \alpha\)」:安定分布を特徴づける数値で,\(\small 0\lt \alpha\leq 2\) の値を取る.
  • 歪度指数「\(\small \beta\)」:分布の歪度を示す数値で, \(\small -1\leq \beta\leq 1\) の値を取る. \(\small \beta=0\) のときは左右対称,\(\small \beta=-1\) のときは右側が切れ, \(\small \beta=1\) のときは左側が切れる形のグラフになる.
  • スケール指数「\(\small \gamma\)」:分布の大きさを決める正の実数値.
  • 移動指数「\(\small \delta\)」:分布の平行移動を決める実数値.
安定分布を \(\small S(\alpha,\beta,\gamma,\delta)\) で表すことにすると, 既知の幾つかの確率分布は,次のような安定分布と対応することが知られている (参照).
  • 正規分布 N\(\small (\mu,\sigma^2)\):安定分布 \(\small S(2,0, \sigma/\sqrt{2},\mu)\) と対応する.
  • コーシー分布\(\small {\rm Cauchy}(x;\gamma,\delta)\) : 安定分布 \(\small S(1,0,\gamma,\delta)\) と対応する.
  • レビ分布 \(\small {\rm Levy}(x;\delta,\gamma)\): 安定分布 \(\small S(1/2,1,\gamma,\gamma+\delta)\) と対応する.
なお,コーシー分布とレビ分布の確率密度関数は次式で定義される. \[\begin{align*}\small {\rm Cauchy}(x;\gamma,\delta)& \small = \frac1{\pi}\frac{\gamma}{(x-\delta)^2+\gamma^2}\\ \small {\rm Levy}(x;\delta,\gamma)& \small=\sqrt{\frac{\gamma}{2\pi}} (x-\delta)^{-\frac32}\exp\left(\frac{\gamma}{2(x-\delta)}\right) \end{align*}\]
▲戻る(トップメニュー)

■確率密度関数
以下では,\(\small \gamma=1, \delta=0\) とする. 歪度指数を \(\small \beta=0\) として, \(\small \alpha=\frac23(黒), 1(赤), \frac32(緑), 2(青)\) の場合の確率密度関数の グラフを描画すると下図のようになる.指数の値が大きくなるにつれ, 頂点部分が下がってくる.シアンは標準正規分布 N(0,1) であり, 安定分布の密度関数では 「dstable(x,2,0,1/sqrt(2),0)」である.

以下は,裾部分 (\(\small x\geq 3\)) を比較したものである. 追加部分のグラフでは,範囲指定は省略できる.

下図では,\(\small \alpha=1\) を固定して,歪度指数を \(\small \beta=-1(黒), -1/2(赤), 0(緑), 1/2(青), 1(シアン)\) に変化させたものである.歪度指数が増えるにつれ,頂点部分が右側から 左側に移動していいく.\(\small \beta=-1\) では裾部分が左側にあり, \(\small \beta=1\) では右側にある.

▲戻る(トップメニュー)

■累積分布関数
下図では,\(\small \beta=0\) の場合の累積分布関数を示した. 横軸は対数軸とした.


▲戻る(トップメニュー)

■stableSlider
上記では個別の関数を指定してグラフを描画したが, 「stableSlider」を利用すると, ボタンをスライドさせて4つのパラメータを インターラクティブに変化させながらグラフを確認することができる. この関数を利用するには,「fBasics」の他に「tclck」も読み込んでおく必要がある.

最初に「stableSlider()」を実行すると,左側にスライドさせるボタンが現れる.

グラフ画面には \(\small \alpha=1.8, \beta=0, \gamma=1, \delta=0\) の場合の確率密度関数と 累積分布関数のグラフが表示される.

ボタンをスライドさせると グラフがインターラクティブに変化する. いろいろと動かしてみることで,各パラメータの意味を把握することができる. 下段の累積分布関数では,下位と上位のそれぞれ1%点の値も表示される. 「pm」ではパラメータの定義の仕方を選択する. 3種類から選択できるが,その意味を把握するのは専門的すぎる. 通常はデフォルトの「pm=0」でよいと思われる. このページの画面はすべて「pm=0」の場合のグラフである.


▲戻る(トップメニュー)

■裾部分の漸近近似
安定分布の裾部分は,ベキ関数で漸近近似することができる. 一般に,安定分布 \(\small S(\alpha,\beta,\gamma,\delta)~(0\lt \alpha\lt 2)\) の 確率密度関数を \(\small f(x)\),累積分布関数を \(\small F(x)\) とすると, \(\small C_{\alpha}=\frac1{\pi}\Gamma(\alpha)\sin\frac{\alpha\pi}{2}\) に対して 次のことが成り立つ(参照, (p.5)). \[\small 1-F(x)\sim (1+\beta)C_{\alpha}x^{-\alpha}\] \[\small f(x)\sim \alpha(1+\beta)C_{\alpha}x^{-(\alpha+1)}\] このことを,実際に確認してみる. 最初に,係数と右辺の関数を定義する.

次に,確率密度関数について確認する.近似の状況を見るため, 横軸は対数軸とした.指数 \(\small \alpha\) の値が大きいほど, 近似の範囲が広いように思える.

\(\small \alpha=\frac23, \beta=0\) のとき.

\(\small \alpha=1, \beta=0\) のとき.

\(\small \alpha=\frac32, \beta=0\) のとき.

このことを逆に見ると,ベキ関数の裾部分を安定分布で近似できることになる. \[\small \frac{\alpha}{x^{\alpha+1}}\sim \frac1{(1+\beta)C_{\alpha}}f(x)\] 下記は,\(\small \alpha=1\) のベキ関数 \(\small 1/x^2\) について, 安定分布 \(\small S(1,\beta,1,0)\) において \(\small \beta=-\frac12(赤), 1(緑), \frac12(青)\) の場合の近似の状況をみたものである. なお,ベキ分布(パレート分布)の確率密度関数を利用するので, 「library(Pareto)」は実行済みとする.


▲戻る(トップメニュー)

■確率変数の和の分布
安定分布の定義によれば, 確率変数 \(\small X\) の独立なコピー \(\small X_1, X_2, \ldots,X_n\) に対して \[\small X_1+X_2+\cdots+X_n\overset{d}{=}a_nX+b_n\] が成り立つような定数 \(\small a_n, b_n\) が存在する. 指数 \(\small \alpha\) の安定分布の場合は \(\small a_n=n^{\frac1{a}}\) である [参照,(8)]. \(\small b_n=0\) のときは厳密に安定であるといわれる. 以下、乱数を利用して,このことを確認する. 和の値は大きくなるので,平均で考えると, \[\begin{align*}\small \frac1{n}(X_1+X_2&\small +\cdots+X_n)\\ &\small \overset{d}{=}\frac{n^{\frac1{a}}}{n}X+\frac{b_n}{n}\end{align*}\] を考えることになる.

最初に,平均からなるベクトルを生成する関数を定義して,たとえば \(\small S(\frac32,0,1,0)\) にしたがう乱数平均からなるベクトルを生成する.

一つにまとめて「summary」を出力させる.

厳密に安定の場合に,右辺の乱数を生成する関数を定義する.

1つにまとめて「summary」を実行する.

刻み幅を 0.2 としてヒストグラム情報を取得する.

下図は,このうち10個平均(x10)のヒストグラムを描画して, 対応する右辺のヒストグラムを折れ線で重ねたものである. ほぼ同じ分布になっていることが分かる.

下図は,\(\small \alpha=\frac23\) の場合である. この場合は裾が広いので,最小値や最大値の値によっては 計算にかなりの時間がかかることがある.


▲戻る(トップメニュー)

■stableFit
パッケージ「fBasics」には,与えられたデータから 安定分布の4つのパラメータを推定する関数が含まれている. 「stableFit」である(参照,(p.30)). 試みに,\(\small S(3/2,0,1,0)\) にしたがう乱数を発生させて, どのようなパラメータとして返されるかを確かめる.

最初に,\(\small [-100,100]\) を 0.2刻みにした値を用意し, その値での確率密度関数の値を「stb0」とする. 次に,要素数1000個の乱数を発生させて「stb1」とする. その値からパラメータがどのように推定されるかを 「stableFit(stb1)」により確かめる. この関数の適用制限は \(\small \alpha\gt 0.6\) になっている.

結果が表示されると同時に,下図のようなグラフも描画される. 与えられたデータから5つの分位数を求めて \(\small \alpha, \beta\) を 評価しているようであるが, パッケージのドキュメントを見ても図の見方については解説されていない.

次に,正しい値であるかどうかをグラフで確認する. 最初に,○印を繋いで確率密度関数のグラフを描き, 次に,「stb1」を0.5刻みとしたヒストグラムを作成する. その上に,推定されたパラメータによる密度関数のグラフを重ねた. それなりのグラフになっていることが分かる.


★ 「fBasics」を読み込むと, 他の分布でもパラメータを推定する関数がある.x をデータとするとき, たとえば,「nFit(x)」で 正規分布の平均と標準偏差の推定とグラフが表示され, 「tFit(x)」でt分布の自由度の推定とグラフが表示される(図は略).
▲戻る(トップメニュー)

■ベキ分布の和の分布
ベキ分布の和」の箇所では, ベキ分布 Pareto(1, 3/2), Pareto(1, 2/3) にしたがう 確率変数の和の分布は安定分布に収束することをみた. 「stableFit」を利用して,その安定分布のパラメータの値を推定してみる.

なお,ここでは,パッケージの「Pareto」と 「gsl」を読み込み済みとし, その箇所で定めた諸関数をそのまま利用することとする. 使用するのは,下記の関数である.

  • 「par_mean」:ベキ分布の平均からなるベクトルを生成する.
  • 「nor_pmean」:そのベクトルを標準化する.
    その定義の中では,係数「coef_a」「coef_b」「coef_C」も使用する.
  • 「f23」「f32」:\(\small \alpha=2/3, 3/2\) の場合の確率密度関数.
    それを定義するために特殊関数「hyperg_U」「w」が使用される.
[1] Pareto(1, 3/2) の場合
下図は,Pareto(1, 3/2) にしたがう乱数を用いて,1000個平均の値からなる 大きさ1000のベクトルに「stableFit」を適用したものである.

下図では,示された値をもとにヒストグラムと安定分布を重ねたものである.

次に,このデータを標準化して同様のことを試みる. 下図は,標準化したデータに「stableFit」を実行したものである. 標準化前の値と比較すると,\(\small \alpha, \beta\) の値は同じであるが, スケール指数と移動指数が異なっている.

標準化したデータをもとに, 階級値の箇所で確率密度関数 \(\small f_{3/2}(x)\) の値を○印で置いて線で繋ぐ. ○印が階級値の箇所に配置されるように,分割する範囲をちょっとずらした. そして,ヒストグラムと推定されたパラメータによる安定分布のグラフを重ねると 次のようになる.

確率密度関数 \(\small f_{3/2}(x)\) の他の3つのパラメータが不明であるが, 以上の推定値とグラフをみると \(\small \alpha=\frac32, \beta=1, \gamma=1, \delta=-1\) であるように思われる. 下図では,\(\small f_{3/2}(x)\) のグラフ上に ○印を置いて線で繋ぎ,dstable(x,3/2,1,1,-1) のグラフを赤で重ねたものである. グラフは一致していることが分かる.


▲戻る(トップメニュー)

[2] Pareto(1, 2/3) の場合
同様にして,\(\small \alpha=\frac23\) の場合の 確率密度関数 \(\small f_{2/3}(x)\) について, 他の3つのパラメータの推定を試みる.

以上の値から,\(\small f_{2/3}(x)\) の他の3つのパラメータは, \(\small \beta=1, \gamma=1, \delta\approx 1.7\) ではないかと思われる. 下図のように,ほぼ重なっている.

▲戻る(トップメニュー)

[3] Pareto(1,1)の場合
このベキ分布の確率密度関数は \(\small 1/x^2\) である. この分布にしたがう確率変数の平均がどのような分布を示すかは, すでに「標本平均の分布」の箇所で調べた. 以下では,標準化した場合のパラメータについて同様の推定を試みる.

以上の値から,Pareto(1,1) にしたがう確率変数の和を標準化した 収束先の安定分布の3つのパラメータは, \(\small \beta=1, \gamma=1, \delta=0\) ではないかと思われる.

標準化する前のデータについても調べると,次のようになる.

★以上の [1]〜[3] の結果をみると,Pareto(1, \(\small \alpha\)) に したがう乱数平均を標準化すると,安定分布 S(\(\small \alpha, 1, 1, \delta\)) に 近づいていくように思われ,\(\small \alpha=1\) のときは \(\small \delta=0\), \(\small \alpha\lt 1\) のときは \(\small \delta\gt 0\),そして \(\small \alpha\gt 1\) のときは \(\small \delta\lt 0\) であるように思われる.


▲戻る(トップメニュー)

レビ分布  [script]

■レビ分布とは
レビ分布は,標準正規分布 N(0,1) にしたがう確率変数を \(\small Z\) とするとき, \(\small 1/Z^2\) のしたがう確率分布として定義され,確率密度関数は次式で表される (参照). \[\small f(x)=\frac1{\sqrt{2\pi}}x^{-\frac32}e^{-\frac1{2x}}\quad (x\gt 0)\] \(\small X\) が N\(\small (\mu,\sigma^2)\) にしたがうときは \(\small Z=(X-\mu)/\sigma\) は N(0,1) にしたがうので, その場合の確率密度関数は次式で表される. \[\small \sqrt{\frac{\sigma}{2\pi}}(x-\mu)^{-\frac32} \exp\left(-\frac{\sigma}{2(x-\mu)}\right)\quad (x\gt \mu)\] この確率分布を Levy(\(\small \mu, \sigma\)) で表すことにすると, この分布は安定分布であり S(\(\small 1/2, 1, \sigma, \mu+\sigma\)) と対応する(参照). したがって,Levy(0,1) は S(\(\small 1/2,1,1,1\)) と対応する. この確率分布は,平均も分散も存在しない.

「R」の標準機能にレビ分布は含まれていないので, この分布を扱うには新たなパッケージを読み込む必要がある. たとえば,「rmutil」(Utilities for Nonlinear Regression and Repeated Measurements Models)というパッケージにレビ分布に関連する関数が登録されている (Document). このパッケージを「パッケージのインストール」を利用してインストールして, 「library(rmutil)」により読み込ませる. また,安定分布のパッケージ(stabledist)と「truehist」を登録している パッケージ「MASS」も読み込んで,以上のことを実際に試してみる.

▲戻る(トップメニュー)

確率密度関数のグラフ
最初に,レビ分布の確率密度関数について確認する. まず,前述の確率密度関数を「levy(x, m, s)」として定義して 「levy(x,0,1)」を点列で描画し,「rmutil」に登録されている レビ分布の確率密度関数「dlevy」のグラフと比較する. 当然ながら,同じであることが確認できる.

次に,同様の点列を,レビ分布に対応する安定分布「dstable(x, 1/2, 1, 1)」 のグラフと確認する.当然ながら,これも一致する. グラフは上図と同じなので省略する.

▲戻る(トップメニュー)

N(0,1)の2乗の逆数の分布
レビ分布は,標準正規分布 N(0,1) にしたがう確率変数 X に対して,1/x^2 の 分布として定義される.このことを,実際に確認してみる. ただし,標準正規分布は 0 の近くの値が多く,その値の2乗の逆数をとると 値が非常に大きくなる.そのため, ヒストグラムの描画には時間がかかり,場合によっては描画できない場合がある. そのため,下記のベクトル「x」は何度か作成して, 最大値が5桁程度の場合で試した方が良い. 下記の値は8回目で生成された.

そのベクトルのヒストグラムを描画して「dlevy(x,0,1)」と重ねると, ほぼ一致していることが分かる.

▲戻る(トップメニュー)

平均の分布
次に,前と同様にして レビ分布が安定分布であることを確認する. レビ分布の指数は \(\small \alpha=\frac12\) であるので, 厳密に安定とすると,\(\small a_n/n=n^{1/\alpha}/n=n\) となるので, \[\begin{align*}\small \frac1{n}(X_1+X_2&\small +\cdots+X_n)\\ &\small \overset{d}{=}nX\end{align*}\] となるはずである. まず,必要な関数を定義する.

次に,5個平均の場合と,対応する右辺の「nX」のベクトルを生成する. 値が大きくなるので最大が6桁程度になるまで何度か繰り返す.

次に,それをもとにヒストグラム情報を取得して, 刻み幅が1の場合の平均の棒グラフと,「nX」の折れ線グラフを重ねると, ほぼ同じ分布になっていることが分かる.

▲戻る(トップメニュー)

copyright