「R」に関する備忘録

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

そのパッケージは自動的に展開されてチェックサムも行われ， 「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」に書いています)．

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

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(緑)\) の場合のグラフです．

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

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

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

下図は，qPareto(x,1,1)のグラフです． 縦軸を対数軸としています．

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

下図では，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」として重ね合わせました．

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

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

1. GNU科学技術計算ライブラリー，リファレンス・マニュアル[578頁]
2. Package 'gsl'

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

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

ベキ分布の確率密度関数は単純ですが， この分布にしたがう確率変数の和の分布は簡単な式で表すことはできません． この分布は指数 \(\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))．

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

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

安定分布は，確率変数 \(\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次式の分布と一致することを示しています． 安定分布は，この関係を平均や分散が存在しないような確率分布に 拡張して考えて得られたものです．

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

• dstable：確率密度関数
• pstable：累積分布関数
• qstable：分位関数
• rstable：乱数生成

• 特性指数「\(\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\)」：分布の平行移動を決める実数値．

• 正規分布 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)\) と対応する．

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

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

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

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

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

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

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

• 「par_mean」：ベキ分布の平均からなるベクトルを生成する．
• 「nor_pmean」：そのベクトルを標準化する．
  その定義の中では，係数「coef_a」「coef_b」「coef_C」も使用する．
• 「f23」「f32」：\(\small \alpha=2/3, 3/2\) の場合の確率密度関数．
  それを定義するために特殊関数「hyperg_U」「w」が使用される．

★以上の [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\) であるように思われます．

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