Maximaを活用した数学学習

in [数ナビの部屋]

「Maxima」を活用した数学学習を取りまとめました.

「Maxima」を活用して数学の学習ロードを駆け抜けよう!

  • MathJaxを使用しているので、 スマホでは表示に時間がかかることがあります。
  • モバイル利用(Android)でのメニュー選択は、 SiteMapを利用するか、 または「長押し」から「新しいタブを開く」を選択してください。
■ 数式処理ソフト「Maxima」を活用した数学学習


[御案内] 「Maxima」(マキシマ)は,フリーの数式処理ソフトです. 有料の Mathematica や Maple に劣らないレベルの数式処理が可能であり, Linux,Windows,MacOSのみならず,Android版もあります. ここでは,数学学習での Maxima の活用法について解説します.

[お知らせ] スマホ(Android)版Maximaの解説本を出版しました. 計算問題やグラフの確認をするときに非常に重宝します. フリーソフトなので一度試してみてください. PC版のコマンドレファランスとしても利用できます。
いつでも・どこでも・スマホで数学!
試し読み 森北出版 amazon 楽天 honto 7net 紀伊國屋電子版読書メーター

■数学学習での活用

以下では,「TeXmacs」+「Maxima」の画面で基本的な使い方を解説します.
詳しい解説は,リンク集に登録したサイトを参照してください.
確率・統計
 統計に特化したソフトウェアとしては,SPSSSASRなどが有名ですが, Maximaにも統計向けのパッケージが備わっています.

■いろいろな確率分布

 いろいろな統計解析を行うには, 特定の統計量がどのような分布をするのかが重要です. Maximaは下記のように多様な分布を扱うことができます. 世の中の多くのことは、 「正規分布」ではなく「ベキ分布」に従っているようです。ベキ分布:リンク集」も作成したので 参照してください。
<Maximaで扱える確率分布>
離散分布
一様分布,二項分布,負の二項分布,ポアソン分布,ベルヌーイ分布,幾何分布,超幾何分布,等々
連続分布
一様分布,正規分布,t 分布,カイ二乗分布, F 分布,指数分布,対数分布,ガンマ分布,ベータ分布, ロジスティック分布,パレート分布,ワイブル分布,レイリー分布,ラプラス分布, コーシー分布,ガンベル分布

 確率分布に関するパッケージ「distrib」を読み込むと, 確率密度関数,分布関数,分布関数の逆関数,平均,分散,標準偏差,歪度係数, 尖度係数,乱数などを個々の確率分布ごとに求めることができます. 詳細は,マニュアルの「52.distrib」を参照してください.


べき分布

 通常の統計の教科書には書かれていないことが多いのですが, 近年,多くの分野で「べき分布」が大きく取り上げられています. これは,確率密度関数が \[\small f(x)=\frac{C}{x^{a}}\quad (x\geq b, b>0)\] の形で表される確率分布で「パレート分布」とも呼ばれます. 一般には,このようなべき乗で表される関係を「べき乗則」といいます. 重力やクーロン力のような関係にあるものもべき乗則にしたがっています. 自然界の現象のみならず, 市町村の人口やWeb上の検索ランキングなど, いろいろな分野でべき乗則にしたがう現象が指摘されており, その背後には何か共通の原理があるのではないかと研究が進められています. 「ベキ分布:リンク集」も参照してください。


確率密度関数
 べき分布の確率密度関数 \(\small f(x)={C}/{x^a}\) の分子の \(\small C\) は規格化定数と呼ばれるものです. この関数が確率密度関数であるためには, \(\small \int_b^{\infty}f(x)\,dx=1\) でなければならないので, \[\begin{align*} \small \int_b^{\infty}f(x)\,dx &\small =\bigg[\frac{C}{(-a+1)x^{a-1}}\bigg]_b^{\infty}\\ &\small =\lim_{x\to\infty}\frac{C}{(-a+1)x^{a-1}}\\ &\small \qquad -\frac{C}{(-a+1)b^{a-1}}\\ &\small =-\frac{C}{(-a+1)b^{a-1}}=1\\ &\small \qquad (a \gt 1)\\ \small \therefore\quad C&\small =(a-1)b^{a-1} \end{align*}\] したがって,\(\small a\gt 1\) のとき,確率密度関数は次のように表せます. \[\small f(x)=\frac{(a-1)b^{a-1}}{x^a}\] 簡単のため \(\small a-1\) を \(\small a\) と置き直して, \[\small f(x)=\frac{ab^a}{x^{a+1}}\] で考えます.ただし,\(\small a>0, b>0, x\geq b\) とします. この式は,次のようにも表せます. \[\small f(x)=\frac{a/b}{(x/b)^{a+1}}\] この関数を, Maximaでは「pdf_pareto(x,a,b)」により参照することができます. 以下では,簡単のため \(\small b=1\) として, \[\small f(x)=\frac{a}{x^{a+1}}\quad (a>0, x\geq 1)\] で考えることにします.「pdf_pareto(x,a,1)」を考えることになります.
 以下では,パッケージ「distrib」と「descriptive」を読み込んだ上で 以上の計算を行っています.

 (%i3)では,無限積分 \(\small \int_{b}^{\infty}c/x^a\,dx\) を 計算しています.\(\small a\) が \(\small -1\) と一致するかどうかが 問われるので「no」を回答します.次に,\(\small a-1\) の符号が問われるので 「pos」と回答すると,結果が表示されます. (%i7)では, \(\small a-1\) から \(\small a\) への置き換えは自分で行っています. この後の計算のため,(%i8)では \(\small a>0\) を仮定します.

 次に,累積分布関数 \(\small F(x)\) を計算してみましょう. \[\begin{align*} \small F(x) &\small =\int_{1}^{x}\frac{a}{t^{a+1}}\,dt\\ &\small =\big[-\frac1{t^a}\bigg]_{1}^{x} =1-\frac1{x^a} \end{align*}\]  このように,べき分布の累積分布関数も べき関数で表されます.この関数は, 「cdf_pareto(x,a,1)」により参照することができます.


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

グラフの比較
 この確率分布のグラフがどのようになるかを, 他の確率分布と比較してみましょう. \(\small x\geq 1\) の部分だけを考えるので, 他の分布も\(\small x=1\) が中心になるように平行移動して, 正規分布は「pdf_normal(x,1,1)」, コーシー分布は「pdf_cauchy(x,1,1)」として比較すると, \(\small a=1\) の場合は次のようになります.

 確率密度関数の形から,正規分布は指数関数的に減少するのに対して, コーシー分布やべき分布はべき関数的に減少していきます. 要するに,減少の仕方が正規分布と比べて緩やかであるということです. 下図では, べき分布 \(\small f(x)={a}/x^{a+1}\) の指数の違いによる グラフを比較しています.

 \(\small a\) の値が大きくなるにつれて減少の度合いが大きくなりますが, そのグラフは \(\small x\) が大きくなると正規分布よりも上側になっていきます.

 参考までに,\(\small b\) の値の違いも見ておきましょう. この値を増やすと,\(\small x\geq b\) であることからグラフは右にずれていきます. そこで,同じ箇所から始まるように平行移動して考えると次のようになり, 値が大きくなるにつれて裾が厚くなるようです.

 累積分布関数は「cdf_pareto(x,a,b)」で表されます. 以下は,\(\small a=1,2,3,~ b=1\) の場合のグラフです.

 以下は,\(\small a=1,~ b=1,2,3\) の場合で, \(\small x=1\) の箇所が中心になるように移動しています.


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

上側%点
 次に,\(\small P(X\geq n)\) の値を求めてみましょう. 累積分布関数 \(\small F(x)\) は「cdf_parete(x,a,b)」で表されるので, この場合は「\(\small 1-{\rm cdf}\_{\rm pareto}(n,a,1)\)」を求めることになります.実際には, \[\small {\rm cdf}\_{\rm parete}(x,a,1)=1-\frac1{x^a}\] であるので, \[\small 1-{\rm cdf}\_{\rm parete}(x,a,1)=\frac1{x^a}\] となります.したがって,\(\small a=1\) のとき, つまり \(\small f(x)=1/x^2\) のときは, \(\small 1-F(x)=1/x\) です. 後でみるように, \(\small a=1\) のとき平均と分散は存在しません.

 上図は,\(\small a=1,2,3\) のそれぞれについて, \(\small P(X\geq n)\) の値を \(\small 2\leq n\leq 10\) の範囲で 求めてグラフ化したものです. 横軸は \(\small n\),縦軸は P(\(\small X\geq n\)) です. 当然ながら,\(\small n\) の値が増えるにつれ \(\small P(X\geq n)\) の 値は小さくなっていきますが,\(\small a=1\) の場合は \(\small P(X\geq 10)=1/10=0.1\) であり,かなり大きな確率です.

 検定では,上側確率や下側確率が重要です.\(\small \alpha\) の 値が与えられたとき, \(\small P(X\geq x)=\alpha\) となる \(\small x\) は どのような値でしょうか.それは \(\small P(X\leq x)=1-\alpha\) でもあるので,累積分布関数の逆関数, つまり「quantile_pareto(x,a,b)」を利用すると \(\small \alpha\) を与えて \(\small x\) の値を求めることができます.

 ここでは小数値を表示する必要があることから,見やすさの観点で, 上図では(%i16)で小数の桁数を4桁に制限しています. (%i17)では,上側確率 \(\small \alpha\) を入力すると, \(\small a=1,2,3\) の場合について 累積確率が \(\small 1-\alpha\) になる値を返す関数を定義しています. (%i18)〜(%i21)では \(\small \alpha=0.5, 0.25, 0.01, 0.005\) の 場合を求めています.実際に検定しようとすると, かなり大きな値まで考える必要があることが分かります.たとえば, \(\small a=1\) の場合の上側1%点は「100」です.

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

対数化変換
 べき分布の確率密度関数 \(\small y=a/x^{a+1}~(a>0, x\geq 1)\) において, 両辺の対数(常用対数)をとると \[\small \log{y}=\log{a}-(a+1)\log{x}\] となります.ここで,2つの座標軸を対数軸にして \(\small Y=\log{y}, X=\log{x}\) とおくと \[\small Y=\log{a}-(a+1)X\] という1次式で表され,\(\small X, Y\) に関するグラフは直線になります. べき分布のグラフは, 急激に減少した後の詳細が把握しにくいので, 対数を取って考えることが多いようです. 確率密度関数を \(\small f(x)=ax^{-(a+1)}\) の形で表したときの 指数部分が直線の傾きになっています. また,両対数グラフが直線 \(\small Y=\log{C}+\beta X\) となる場合は, \[\small \log{y}=\log{C}+\log{x^{\beta}} =\log{Cx^{\beta}}\] となることから, \(\small y=Cx^{\beta}\) というベキ関数になることも分かります.
 以下は,\(\small a=1,2,3\) の場合のグラフを両対数グラフで描画したものです. \(\small a\) の値が大きくなるにつれ下降が急になることが分かります.

 正規分布(pdf_normal(x,1,1))やコーシー分布(pdf_cauchy(x,1,1))とも 比較してみましょう.下図は,ベキ分布(pdf_pareto(x,1,1))と区間 \(\small [1,5]\) での両対数グラフでの比較です. 正規分布は急激に減少してしまうのに比べて, ベキ分布やコーシー分布は右下がりの直線的に減少しています.

 ここで,べき分布(パレート分布)にしたがう乱数を発生させてみましょう.

 (%i9)で,\(\small a=1, b=1\) のときのパレート分布にしたがう乱数 \(\small n\) 個からなるリストを作成するコマンドを定義します. (%i10)では,要素数1000個のリストを作成して「pare」とし, (%i11)〜(%i13)で大きさ・最小値・最大値を確認しています. 「smin」「smax」はリスト内の要素の最小値と最大値を表示します. 引数はリストか行列で,行列の場合は各列の最小値・最大値を返します. 「lmin」「lmax」も同様の機能を持ちますが, 引数はリストか集合です.

 (%i14)では,生成したリスト「pare」から, 区間 \(\small [1,1001]\) を4000等分したときの 度数分布を得ています. 最大値の値から実際には[1,745]で考えれば十分ですが, ここでは値の切れの良さで1000としています. 「continuous_freq」では区分点[1]と度数分布[2]の2つのリストが返されるので, [2]を指定することで度数分布のリストを得ています. なお,「TeXmacs+Maxima」では(%i13)で最大値が1001を超えるものがあると エラーになるので, 範囲内に納まることを確認してから度数分布を作成しました. このエラーは,wxMaximaでは修正されています.
 (%i15)では,度数分布をそのまま離散グラフとして表示させるため, \(\small x\) 軸側の分点リストを作成しています. (%i16)では,離散グラフとパレート分布の確率密度関数を表示させています. 縦のスケールを揃えるため,密度関数には分割幅と総度数を掛けています.  座標軸を対数軸とするには,「logx」「logy」を指定するだけです.


 このように,急激に減少するものの,そのままダラダラと度数が途切れなく 続く分布が「べき分布」(パレート分布)です. 両座標軸を対数軸に取ると直線的になるのが特徴です. このような特徴をもつ分布は身の周りに多数散見されることから, 世の中で重要な分布は「正規分布ではなくてべき分布なのではないか?」 とも言われているようです.後で述べる 「80:20の法則」や 「ジップの法則」も, 要するに「べき分布」です.


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

平均と分散
 べき分布の平均と分散を計算してみましょう. 確率密度関数を \(\small f(x)=a/x^{a+1}~(a>0, x\geq 1)\) とすると, 次のように計算されます. \[\begin{align*} \small \int_{1}^{\infty}xf(x)\,dx &\small =\int_{1}^{\infty}\frac{ax}{x^{a+1}}\,dx\\ &\small =\int_{1}^{\infty}ax^{-a}\,dx\\ &\small =\big[\frac1{1-a}ax^{1-a}\big]_1^{\infty}\\ &\small =\lim_{x\to\infty}\frac{a}{(1-a)x^{a-1}}\\ &\small \quad -\frac{a}{1-a}\\ &\small =\frac{a}{a-1}\quad (a\gt 1 のとき)\\ \small \therefore\quad E(X)&=\small \frac{a}{a-1}\quad (a\gt 1) \end{align*}\]  上の計算からも分かるように,この積分は \(\small 0\lt a\leq 1\) のときは発散してしまいます. \(\small a\gt 1\) であれば平均が存在します. \[\begin{align*}\small \small \int_{1}^{\infty}x^2f(x)\,dx &\small =\int_{1}^{\infty}\frac{ax^2}{x^{a+1}}\,dx\\ &\small =\int_{1}^{\infty}ax^{1-a}\,dx\\ &\small =\big[\frac{a}{2-a}x^{2-a}\big]_1^{\infty}\\ &\small =\lim_{x\to\infty}\frac{a}{(2-a)x^{a-2}}\\ &\small \quad -\frac{a}{2-a}\\ &\small =\frac{a}{a-2}\quad (a>2 のとき)\\ \small \therefore\quad V(X)&\small =E(X^2)-E(X)^2\\ &\small =\frac{a}{a-2}-\left(\frac{a}{a-1}\right)^2\\ &\small =\frac{a}{(a-1)^2(a-2)} \end{align*}\]  ここでも,\(\small a\leq 2\) であれば分散は存在しません. 同様に計算していくと,\(\small a\leq 3\) のときは歪度が, \(\small a\leq 4\) のときは尖度が存在しません.

 次に,「分散が存在しない」ということを, コーシー分布のときと同様の手法で確認してみましょう. 標本分散だと値が大きくなるので,標本標準偏差で確認します.

 (%i19)では,乱数を発生させる個数を横軸にするため, その分点を指定するリストを作成しています. 10刻みで1000個まで考えることにします. (%i20)では,\(\small a=1.5\) の場合のべき分布(パレート分布)にしたがう乱数が \(\small n\) 個からなるリストを作成するコマンドを定義しています. (%i21)(%i22)は,そのリストの標本平均(mean)と標本標準偏差(std)からなるリストです. 計算する個数は, 10個からはじめて10刻みに増やしながら1000個まで計算しています.

 \(\small a=1.5\) のときの乱数なので,  この確率分布の平均は存在しても分散は存在しません. 平均は \(\small \frac{a}{a-1}=\frac{1.5}{0.5}=3\) です. グラフを見ると,標本平均は(この場合は)ずっと一定の値を保っていますが, 標本標準偏差の値は個数を増やしても安定しません. 乱数なのでグラフを描画させるごとに異なるグラフが描画されます. 標本平均のグラフと標本標準偏差のグラフでは異なる乱数が用いられています. 標本標準偏差が100を超える場合もありました。
 以下は,標本平均と標本標準偏差の変化を個別にみたものです. 上段は標本平均の変化,下段は標本標準偏差の変化です. いずれも,上図とは異なる乱数によるものです. 

 これをみると,個数を増やすと(確率分布としての平均が存在するので) 標本平均の値の変動は小さくなっていきますが, 標本標準偏差の値の変動は逆に大きくなっていくようです.
 標本標準偏差が大きな値を取るとき, 具体的にはどのような値になっているのでしょうか. 試みに「plist2(10)」により10個の成分を持つリストを発生させ, 標本標準偏差が大きな値になるまで何度か試したところ, 次の例が生成されました.

 実際には,10回以上試してこのリストが得られました. 整数部分だけを書き出すと, 「\(\small 1, 6, 493, 2, 1, 1, 25, 11, 1, 3\)」 という値でした.「493」が一つ混じったことが原因のようです.

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

80:20の法則
 パレート分布は,イタリアのヴィルフレド・パレートにより提唱されたものです. 経済社会における所得分布を分析して, 2割の高額所得者のもとに社会全体の8割の富が集中して, 残りの2割の富が8割の低所得者に配分されるという「パレートの法則」を 見いだしました.この法則は「80:20の法則」と呼ばれることもあり, 経済学のみならず自然現象や社会現象などの様々な事例にあてはまるとされます. Wikipediaでは, 下記の例が紹介されています.
  • ビジネスにおいて、売上の8割は全顧客の2割が生み出している。よって売上を伸ばすには顧客全員を対象としたサービスを行うよりも、2割の顧客に的を絞ったサービスを行う方が効率的である。
  • 商品の売上の8割は、全商品銘柄のうちの2割で生み出している。
  • 売上の8割は、全従業員のうちの2割で生み出している。
  • 仕事の成果の8割は、費やした時間全体のうちの2割の時間で生み出している。
  • 故障の8割は、全部品のうち2割に原因がある。
  • 住民税の8割は、全住民のうち2割の富裕層が担っている。
  • プログラムの処理にかかる時間の80%はコード全体の20%の部分が占める。
  • 全体の20%が優れた設計ならば実用上80%の状況で優れた能力を発揮する。
 他にも,あります.リスク管理・Navi では次の事象があげられています. ベキ分布の実例も参照してください.
  • ガラスを床に落とした際の破片の大きさの分布
  • 地震の大きさと発生頻度
  • 山火事の被害面積ごとの発生頻度
  • 株価や為替などの市場価格の変動
  • 所得や純資産の分布
  • 本の売上げ分布
  • 論文の数と引用された回数の関係
  • 戦争の発生頻度と死者数
 そこで,べき分布 \(\small f(x)=a/x^{a+1}~(x\geq 1)\) において, \(\small x\) を富とした場合を, 「須鎗弘樹:複雑系のための基礎数理(牧野書店),p9」を もとに考えてみましょう.
 富が \(\small x\) 以上である人の割合 \(\small S\) は 既に計算済みの式から \[\small S=1-P(X\lt x)=\frac{1}{x^a}\] であり,その人達の富は, \[\small \int_x^{\infty} tf(t)\,dt=\frac{a}{(a-1)x^{a-1}}\] で表されます.一方,富全体は,\(\small a>1\) の場合は \[\small \int_1^{\infty} tf(t)\,dt=\frac{a}{a-1}\] となるので,富全体の中で,富が \(\small x\) 以上の人達の富の 占める割合 \(\small T\) は \(\small T=1/x^{a-1}\) で表されます. つまり,富が \(\small x\) 以上である人達の割合 \(\small S\) と, その人達の富が全体の富の中で占める割合 \(\small T\) の間には \[\small T=S^{\frac{a-1}{a}}\] という関係があることになります. したがって,富の80%を20%の人達が占めているという関係は, \[\small T:S=80:20\] ということから \[\small S^{\frac{a-1}{a}}:S=80:20\]  これより \(\small 4=S^{-\frac1{a}}\) が得られ,両辺の対数をとることにより \[\small a=-\frac{\log{S}}{\log{4}}=-\frac12\log_2{S}\] という関係が得られます. 下図は,\(\small S=P(X\geq x)\) の値と指数 \(\small a\) との関係を 図示したものです.

 指数 \(\small a\) の値は,実際には \(\small 1\sim 2\) の値を取ることが 多いようなので,\(\small S=P(X\geq x)\) の値は \(\small 0.05\sim 0.25\) の 範囲にあります.つまり,全体の5%から25%の人達により, 全ての富の80%が占められているということになります.

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

ジップの法則
 ベキ分布の特別な場合に「ジップ(Zipf)の法則」と呼ばれるものがあります. これは,いろいろな現象において, 「ある事象 \(\small A\) の規模と事象 \(\small A\) の順位の積は一定である」 という経験則です. George K. Zipfは言語学者で,小説「ユリシーズ」で使われた 260,430単語に 現れる29,899単語の出現頻度を調べ, 単語の出現回数と出現順位の積が一定であることを見いだしました. つまり,ある単語が \(\small N\) 回現れ, その単語の出現順位が \(\small n\) 番目とすると \(\small Nn=C~(定数)\) であること, つまり,\(\small N={C}/{n}\) という関係があることを見いだしました.
 ある事象の総量と順位との関係では,他にも多数のものが ジップの法則にしたがうとされています. Wikipediaでは, 次のような例が挙げられています. ベキ分布の実例も参照してください.
  • 日本の都市の人口
  • 地震の規模
  • ホームページのリンク数
  • 企業の所得額
  • 音楽における音符の使用頻度
  • 細胞内の遺伝子の発現量
  • 固体がわれたときの破片の大きさ

 この法則は,一見するとべき分布 \(\small {C}/{x^a}\) の \(\small a=1\) の場合のように見えるのですが, 以下にみるように,べき分布 \(\small {C}/{x^2}\) の 累積分布関数として捉えることができるようです.
 最初に, ランキングサイトをもとに, 平成27年の国勢調査による県庁所在地47都市の上位20位までの人口と 順位のグラフを以下に示します. 下図は,上位20都市の人口と順位のグラフ,そして両対数グラフです. 縦軸の人口の単位は「千人」としています.

 次に,[丸田]をもとに, 人口 \(\small x\) と順位 \(\small R\) の関係をみてみましょう. まず,上図にあるように, 両対数グラフで考えると,人口と順位との間には直線の関係が見られるます. したがって,対数化変換の箇所で見たように, この2つの変量の間には \[\small x=aR^{-\beta}~(a:正の定数)\] という関係があります. この式より, \[\small R=a^{\frac1{\beta}}x^{-\frac1{\beta}}\] が得られます. ここで,都市総数を \(\small N\) として,\(\small r=R/N\) とおくと, \(\small r\) は全都市の中での順位比率を表します.このとき, \(\small R=rN\) より,上式は \[r=\frac{A}{N}\cdot x^{-\frac1{\beta}}\quad (A=a^{\frac1{\beta}})\] となります.\(\small r\) の定め方から \(\small 0\leq r\leq 1\) であり,この値が小さいほど人口が多いことを, 値が 1 に近いほど人口が少ないことを表します. つまり, 都市の人口を表す確率変数を \(\small X\) とすると, \(\small r\) の値は \[\small P(X\gt x)=r\] を意味する値として捉えることができます. したがって,その逆を考えると, \[\begin{align*} \small P(X\leq x) &\small =1-r\\ &\small =1-\frac{A}{N}\cdot x^{-\frac1{\beta}} \end{align*}\] となります.これは,累積分布関数です. この関数を \(\small F(x)\) とおいて微分することにより, 確率密度関数 \(\small f(x)\) は \[\begin{align*} \small f(x) &\small =F'(x)\\ &\small =\frac{A}{N\beta}\cdot x^{-\frac1{\beta}-1}\\ &\small =\frac{A/(N\beta)}{x^{\frac1{\beta}+1}} \end{align*}\] となります.\(\small \beta=1\) の場合がジップの法則です.


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

分布の相似性
 べき分布は「自己相似性」を持つとされます. それは,スケール変換をしても同じべき乗になることを指します. 実際,べき分布の確率密度関数を \(\small f(x)\) とすると, \[\small f(x)=\frac{a}{x^{a+1}}\] です.今,\(\small x\) を \(\small c\) 倍すると, \[\small f(cx)=\frac{a}{(cx)^{a+1}} =\frac{a/c^{a+1}}{x^{a+1}}\] となり,同じベキを持ちます. この式をグラフでみると,横軸を \(\small c\) 倍したときは, 縦軸を \(\small c^{a+1}\) 倍すると同じグラフになることを示します.

 また,このような性質を持つ関数はベキ関数に限られることを導くことが できます.以下は,「須鎗弘樹:複雑系のための基礎数理(牧野書店),p10」を もとにしたものです.
 \(\small f(x)\) を微分可能な正値関数とします.そして, 任意の \(\small c\) に対して関数 \(\small g(x)\) が存在して \[\small f(cx)=g(c)f(x)\] が成り立ったとします.\(\small x=1\) のときは \(\small f(c)=g(c)f(1)\) となるので,上式は \(\small f(x)\) だけを用いて \[\small f(cx)=\frac{f(c)}{f(1)}f(x)\] と表すことができます.両辺を \(\small c\) で微分すると \[\small xf'(cx)=\frac{f'(c)}{f(1)}f(x)\] ここで,\(\small c=1\) のときを考えると \[\small xf'(x)=\frac{f'(1)}{f(1)}f(x)\] これを変形すると \[\small \frac{f'(x)}{f(x)}=\frac{f'(1)}{f(1)}\frac1{x}\] となります.これを積分すると, \[\small \log{f(x)}=\gamma \log{x}+C,\quad \gamma=\frac{f'(1)}{f(1)}\] これを整理することで,次のようなベキ関数が得られます. \[\small f(x)=e^{C}x^{\gamma}\]

 次に,この相似性をグラフで確認してみましょう. 最初に,ベキ分布(pdf_pareto(x,1,1))の区間 \(\small [3,4]\) の箇所を 単位ステップ関数を利用して切り出して \(\small f(x)\) とします.

 ここでは,\(\small a=1\) の場合なので,たとえば \(\small c=2\) とすると \(\small c^{a+1}=2^2=4\) となります. 関数 \(\small f(x)\) の \(\small x\) を 2 倍して \(\small y\) を 4 倍した グラフを描画すると下図のようになり,もとのグラフとすっかり重なってしまいます. \(\small x\) を2倍すると, グラフの上では \(\small x\) 軸方向に 1/2 倍されます.

 下図は \(\small c=1/2\) の場合です.\(\small [5,10]\) の範囲を示しました. グラフは,\(\small x\) 軸方向に2倍に引き伸ばされます.


 この自己相似性を別な角度から見てみましょう. \(\small a=1\) の場合に「80:20の法則」を考えます. p分位数のコマンド「quantile_pareto(p,1,1)」を利用すると, 下図にあるように [1,5] の範囲に全体の 4/5 が含まれることが分かります. では,残りの \(\small [5,\infty)\) のうちの 80% の箇所はどこでしょうか. それは,全体では \(\small 0.8+0.2\cdot 0.8=0.96\) の箇所なので,同じコマンドを利用すると [5,25] の範囲です. 同様にすると,残りの \(\small [25,\infty)\) のうちの 80%の箇所は,全体では \(\small 0.96+0.04\cdot 0.8=0.992\) の箇所なので [25,125], その残りの \(\small [125,\infty)\) のうちの 80% を占めるのは, 全体では \(\small 0.992+0.008\cdot 0.8=0.9984\) の箇所なので [125,625] の範囲です.範囲が,いずれも 5 のべき乗になっています.

 以上のことを確率密度関数 \(\small f(x)\) を積分する形で考えると, 次のようになります.\(\small a=1\) の場合を考えているので \[\begin{align*} \small \int_{\alpha}^{\beta}f(x)\,dx &\small =\int_{\alpha}^{\beta}\frac1{x^2}\,dx\\ &\small =\bigg[-\frac1{x}\bigg]_{\alpha}^{\beta} =\frac1{\alpha}-\frac1{\beta} \end{align*}\] \[\small \int_1^5f(x)\,dx=\frac45=0.8\] \[\small \int_5^{25}f(x)\,dx=\frac{4}{25}=0.16\] \[\small \int_{25}^{125}f(x)\,dx=\frac{4}{125}=0.032\] \[\small \int_{125}^{625}f(x)\,dx=\frac{4}{625}=0.0064\]

 ここで,2番目の積分式 \[\small \int_5^{25}f(x)\,dx\] で \(\small x/5=t\) という変数変換を行うと, \[\small \int_5^{25}f(x)\,dx=\int_1^5 f(5t)\cdot 5dt\] となります.他の場合も同様であり, 残りの箇所で80%を占める範囲は, スケール変換するといずれも最初の積分範囲と同じ範囲で表されることになります.
 同じことを \(\small f(x)=a/x^{a+1}\) で考え, \[\small \int_1^{c}f(x)\,dx=p\] とします.\(\small p\) の値は,実際には \[\begin{align*} \small p &\small =\int_1^{c}\frac{a}{x^{a+1}}\,dx\\ &\small =\bigg[-\frac1{x^a}\bigg]_1^c\\ &\small =1-\frac1{c^a}\\ \small \therefore\quad \frac1{c^a} &\small =1-p \end{align*}\]  ここで, 残りの \(\small [c,\infty)\) の箇所で同じ割合 \(\small p\) を占める部分は 全体では \(\small (1-p)p=p/c^a\) です. この値は,\(\small x/c=t\) と変数変換して 次の積分を計算しても得られます. \[\begin{align*} \small \int_1^{c}f(ct)\,cdt &\small =\int_1^c\frac{a}{(ct)^{a+1}}\,cdt\\ &\small =\frac1{c^a}\int_1^c\frac{a}{t^{a+1}}\,dt\\ &\small =\frac{p}{c^a} \end{align*}\]  他の箇所でも同様であり, スケール変換すると同じ積分範囲で計算することができます. このように,部分が全体と同じような性質を持つことから, この性質は「スケールフリー」や「フラクタル性」と呼ばれます.

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

平均の分布
 べき分布では,指数の値によっては平均や分散を持たない場合があります. 中心極限定理は 平均や分散が存在するような確率分布に関する定理です. それらが存在しない確率分布では,標本平均はどのような分布になるのでしょうか.
 ここで,ベキ分布で \(\small a=1, b=1\) の場合の 乱数「random_pareto(1,1)」を利用して, 標本平均の分布がどのようになるかを確かめてみましょう. この分布は,平均も分散も存在しません. 平均の存在しない確率分布で乱数を発生させると, コーシー分布の箇所でみたように 標本数を幾ら増やしても標本平均の値は安定しません. その振れ幅は非常に大きなものになります. Maximaのバージョンの古い「TeXmacs+Maxima」では 度数分布の作成で支障をきたすので, 以下ではwxMaximaを利用します. 行番号はリセットして「distrib」と「decriptive」のパッケージを 読み込んでいます.

 (%i3)では,\(\small a=1, b=1\) のパレート分布にしたがう乱数 \(\small n\) 個からなるリストを作成するコマンドを定義しています. (%i4)で,そのリストの平均からなる \(\small N\) 個のリストを作成します. (%i5)(%i6)では, それぞれ2個と5個の標本平均からなる要素数10000個のリストを作成しています. (%i7)(%i8)では,区間 [1, 21] の範囲での 度数分布を作成しています. (実際には [1, 最大値] の範囲で考える必要がありますが, ベキ分布にしたがう乱数を 10000個も発生させると, その最大値はとてつもなく大きな値になるので, ここでは [1, 21] で考えることにします.) 100区分にしたので,0.2刻みの度数分布になります. (%i9)では,グラフを描く際の横軸の分点として, 階級の中点のリストを作成しています., (%i10)(%i11)は,plot2dで描画するときの 離散グラフの書式を定めています.ヒストグラムになるように, 階級幅と指定した区間での総計で割っています.
 以上の準備のもとで, グラフは(%i12)により下図のようになります.

 以下は,同じグラフを両対数グラフで描画したものです.

 これをみると, 平均を取る個数を増やしても正規分布に近づく気配は全くありません. しかし,ある箇所から先では もとの分布と同じような分布の仕方であることに気づきます. 両対数グラフをみると最初の部分を除くと直線状になるので, ある箇所から先ではベキ分布になっていると思われます. 平均を取る個数を増やして試してみても同様であり, 一つの盛り上がりの後を経た後はダラダラとベキ分布的に減少していきます. \(\small n\to\infty\) のときは「安定分布」 と呼ばれる確率分布に近づいていくことが知られています.
 下記は,同様にして乱数の10個平均のリスト「pare10」を作成して固定したとき, その生データをみたものです.横軸は何回目の標本平均であるかを表します. 上から順に,10個平均の1万個のデータ,そのうちの4001〜5000番目のデータ, そして4601〜4700番目のデータです. 縦のサイズは調整しましたが,どの場面でみても 同じような上下動が繰り返されています. 分布としての平均や分散が存在しないので, 標本平均を何回取ってもどのような値になりそうなのか「全く分からない」ということです. おまけに,1万個のデータを見れば分かるように, 10個平均の値が1000を超えることが何度かあります. ときどき,「とんでもないことが起きることがある」ということです.


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

標本平均の差の変化
 前項で,標本平均を求める度に値が大きく変動している様子を見ました. その変動幅はどの程度のものなのかを調べてみましょう.

 後からやり直しているので入力番号がずれていますが, (%i10)では10個平均の値が10001個からなるリストを作成しています. そして,(%i11)では隣同士の差の絶対値からなるリストを作成しています. 以下では,縦軸を対数軸にして,その絶対値がどのように変化しているかを みたものです.激しく振動しているのが分かります.

 次に,この振動の度数分布を調べてみましょう.

 (%i14)では,区間 \(\small [0,100]\) を 500等分したときの 度数分布を得ています.(%i15)は各階級の中央値のリストです. (%i16)でグラフの諸式を定めます.
 以上のもとでグラフを描くと次のようになります.

 変動の絶対値は,ベキ分布になっているようです. 下図は両対数グラフにしたものです. 出だし部分を除くと直線状になっているのが分かります.


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

標本平均の比の変化
 前項で,10個平均の差の絶対値がべき分布になっているのを見ました. ここでは,隣同士の比がどのように変化するかを眺めてみましょう. 最初に,(%i20)により10個平均のリストを作り直します.

 隣同士の比を取る関係で,10001個のリストになっています. (%i22)で,1つ前の標本平均との比で構成されるリストを「pareRatio」とします. 最小値は \(\small 6.48\times 10^{-4}\), 最大値は \(\small 2165.4\) です.

 (%i26)(%i27)で離散グラフの書式を定めておいてから(%i28)でグラフを 描画しています.値が大きすぎるので,縦軸を対数軸(常用対数)に指定しました.
 下図は,一つ前の標本平均との比のグラフです.これも対数軸にしています. 激しく振動していることが分かります. 

 次に,この振動の度数分布を作成してヒストグラムを描いてみましょう.

 (%i31)では,\(\small \log({\rm pareRatio}[i+1]/{\rm pareRatio}[i])\) の 値からなる1万個のデータの度数分布を作成しています. なお,Maximaの「log」は自然対数です.比の値の最小値と最大値は (%i23)(%i24)で確認済みなので, 全ての値が \(\small [-10,10] \) の中に納まります. 実際には \(\small [-8,8]\) で十分ですが,区切りの良い区間にしました. その区間を100等分したので0.2刻みの度数分布です. 以上のもとでグラフを描画したのが下図です.

 参考までに標準正規分布(赤色)も描画しましたが, 平均の比の対数(自然対数)は,なぜかラプラス分布になるようです. ラプラス分布の確率密度関数は,Maximaでは次の式で定義されています. \[\small {\rm pdf\_laplace}(x,a,b)=\frac1{2b}e^{-\frac{|x-a|}{b}}\]  なお,株価の変動でも,価格の比の対数はラプラス分布になるようです [参照(p6)] .

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

分散の有無による比較
 今度は,分散の存在の有無により 標本平均の分布がどのように異なるかをみてみましょう. 「平均と分散」の箇所でみたように, ベキ分布 \[\small f(x)=\frac{a}{x^{a+1}}\] は,\(\small a=1\) のときは平均と分散が存在しません. \(\small a=2\) のとき平均は \(\small E(X)=2\) ですが 分散は存在しません.\(\small a=3\) のときは 平均は \(\small E(X)=\frac32\) , 分散は \(\small V(X)=\frac34\) になります.

 ここでは,行番号をリセットしてやり直しています. (%i3)(%i4)で \(\small a\) の値も変えられるようにして, (%i5)〜(%i7)で10個平均の値が1万個からなるリストを生成しています. (%i9)〜(%i11)でグラフの書式を定めます.
 以上のもとで, \(\small a=1\) の場合の標本平均の分布の様子を描画したのが下図です.

 これは,単に \(\small n\) 回目の標本平均の値を図示したものであり, 縦軸は10個平均の値です.横軸は「回目」を表します. 縦軸を対数軸にとっています. 分布としての平均も分散も存在しないので, 10個平均の値は大きく変動しています. 標本平均が10000を超える場合もあります.
 下図は,\(\small a=2\) の場合です.

 この場合は,平均が 2 で分散は存在しません. 分散が存在しなくても平均が存在することから, 値の変動幅が大幅に狭まっています.
 下図は,\(\small a=3\) の場合です.

 この場合は,平均が 1.5 で分散は 0.75 です. 値の変動幅がさらに狭まっているのが分かります.

 今度は,以上のリストのヒストグラムを見てみましょう.

 (%i15)〜(%i17)では,区間 \(\small [1,21]\) を100等分したときの 度数分布を求めています. 0.2刻みになるので,(%i18)で横軸の分点リストを作成し, (%i19)〜(%i22)でグラフ描画のための書式を設定します.
 以上のもとで,この区間のヒストグラムは下図のようになります. 両対数軸のグラフも示しました.

 平均が存在すると,その箇所での盛り上がりができ, 分散が存在すると平均の近くへの集中が高まっていることが分かります.

 以上の例から,平均や分散が存在しないような確率分布では 標本平均は大きく変動するものの, そのヒストグラムは正規分布とは異なる何か別の確率分布に したがうのではないかと考えられます. それが,別ページで述べる「安定分布」です. 正規分布もコーシー分布も安定分布です. そして,ベキ分布 \(\small a/x^{a+1}\) にしたがう確率変数の和の分布は, 標本数を大きくしていくと,同じ指数 \(\small a\) を持つ 安定分布に収束していくことが証明されています (参照,p.12(19)). そのことを確かめるには安定分布を理解する必要がありますが, 安定分布の確率密度関数を初等関数で表すことができるのは, 正規分布,コーシー分布,そして次の項目にあるレビ分布に限られます. 他の場合は数値的に調べるしかありません. Maximaには「安定分布」に関するパッケージはないので, 以上のことを確かめるため, 安定分布に関するパッケージを持つ「R」を使用することにします. この件に関する詳細は「こちら」を参照してください.


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

べき分布の実例
 ここでは,ベキ分布・自己相似性・べき乗則を持つ 実例について詳細に述べている記事や論文を集めてみました. すでに「パレートの法則」や「ジップの法則」の箇所で 紹介した内容もありますが,具体的なグラフや解説が加えられており, 両対数軸で考えると右下がりか右上がりの直線上に乗っています.


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

レビ分布

確率密度関数
 レビ分布は, \(\small Z\) をN(0,1)にしたがう確率変数とするとき, \(\small 1/Z^2\) のしたがう確率分布です. なぜならば,\(\small U=1/Z^2\) とおいて, \(\small U\) の累積分布関数を考えます. \[\begin{align*} \small P(U\leq u) &\small =P\left(\frac1{Z^2}\leq u\right)\\ &\small =P\left(Z^2\geq \frac1{u}\right)\\ &\small =P\left(Z\leq -\frac1{\sqrt{u}}\right)\\ &\small \qquad +P\left(Z\geq \frac1{\sqrt{u}}\right)\\ &\small =2P\left(Z\geq \frac1{\sqrt{u}}\right)\\ &\small =2\left\{1-P\left(Z\leq \frac1{\sqrt{u}}\right)\right\} \end{align*}\] したがって,N(0,1)の確率密度関数を \[\small \varphi(x)=\frac1{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\] として,\(\small U=1/Z^2\) のしたがう確率分布の 確率密度関数を \(\small f(u)\) とすると \[\small \int_{-\infty}^{u}f(t)\,dt =2\left(1-\int_{-\infty}^{\frac1{\sqrt{u}}}\varphi(t)\,dt\right)\] となります.両辺を \(\small u\) で微分すると, 微積分の基本定理により \[\begin{align*} \small f(u) &\small =-2\frac{d}{du}\left(\int_{-\infty}^{\frac1{\sqrt{u}}} \varphi(t)\,dt\right)\\ &\small =-2\varphi\left(\frac1{\sqrt{u}}\right) \cdot\left(\frac1{\sqrt{u}}\right)'\\ &\small =-2\cdot\frac1{\sqrt{2\pi}}e^{-\frac1{2u}}\cdot \left(-\frac1{2u^{\frac32}}\right)\\ &\small =\frac1{\sqrt{2\pi}}\frac1{u^{\frac32}}e^{-\frac1{2u}}\\ \small \therefore\quad &\small f(x)=\frac1{\sqrt{2\pi}}\frac1{x^{\frac32}}e^{-\frac1{2x}} \end{align*} \] \(\small X\) が N\(\small(\mu,\sigma^2)\) にしたがうときは, \(\small Z=(X-\mu)/\sigma\) とおくと \(\small Z\) はN(0,1)にしたがうので, 同様にすると次の式が得られます. \[\small \sqrt{\frac{\sigma}{2\pi}}\frac1{(x-\mu)^{\frac32}} \exp\left(-\frac{\sigma}{2(x-\mu)}\right)\]

諸関数の定義
 Maximaにレビ分布は登録されていないので, 最初に必要な関数を定義しておきます. 都合上,記号を \(\small \gamma,\delta\) に変更すると, レビ分布の確率密度関数の一般形は \[\small \sqrt{\frac{\gamma}{2\pi}} \frac1{(x-\delta)^{\frac32}}\exp\left(-\frac{\gamma}{2(x-\delta)}\right)\] です.これを Levy(\(\small \delta, \gamma\)) と表すことにします. \(\small \delta\) は平行移動に関わるパラメーターです. 最初に,\(\small \gamma\) の違いでグラフがどのように変わるかを 見てみましょう.wxMaximaを使用します.
 なお,以下のコマンドでは「levi」としました. 後で「levy」とすべきことに気づきましたが, 修正箇所多数なのでそのままにします.

 行番号を初期化して,2つのパッケージを読み込んでおきます. (%i4)でグラフを比較すると,\(\small c\) の値が大きくなるにつれ, 出だしの部分こそ大きく異なりますが, その後の裾の状況にあまり大きな違いはないように見えます. 裾部分の違いを, 「グラフの比較」の(%i9)〜(%i26)の 箇所と同様にして調べてみましょう. \(\small P(X\geq x)\) の値を, 「cdf_levi(x,c):=ldefint(f(t,c),t,0,x)」として 「1\(\small -\)cdf_levi(x,c)」により求めて, 折れ線でグラフ化すると次のようになります.

 これを見ると,確率密度関数の上では大きな違いは見えませんでしたが, 上位確率の値をみると大きな違いが見られます.特に, \(\small c\geq 2\) の場合は,半分以上が 3 以上の裾部分にあります.

 次に,簡単のため \(\small \gamma=1, \delta=0\) の場合を考え, Maximaの書式に合わせて \[\small {\rm pdf\_levi}(x)=\frac{1}{\sqrt{2\pi}} x^{-{\frac32}}e^{-\frac{1}{2x}}\] が確率密度関数の場合を考えます. この分布にしたがう乱数を発生させるには, 累積分布関数の逆関数 を考える必要があります. レビ分布の累積分布関数は, 下記により計算することができます. \[\begin{align*} \small {\rm cdf\_}&\small {\rm levi}(x)=P(X\leq x)\\ &\small =\int_{0}^{x}{\rm pdf\_levi}(t)\,dt\\ &\small =\frac1{\sqrt{2\pi}} \int_{0}^{x}t^{-{\frac32}}e^{-\frac{1}{2t}}\,dt\\ &\small =1-\frac1{\sqrt{2\pi}} \int_{x}^{\infty}t^{-{\frac32}}e^{-\frac{1}{2t}}\,dt\\ \end{align*}\]  しかしながら,この積分は初等積分可能ではありません. 右辺の第2項は,フーリエ解析の箇所で現れた 第ニ種不完全ガンマ関数 に似ています.第ニ種不完全ガンマ関数は, \[\small \Gamma(a,x)=\int_{x}^{\infty}t^{a-1}e^{-t}\,dt\] で定義される関数です.
 ここで,形を揃えるために \(\small 1/2t=s\) として置換積分を行うと, \(\small t=1/2s,~dt=-1/2s^2\) であることから, 次のように変形されます. \[\begin{align*} \small {\rm cdf\_}&\small {\rm levi}(x)\\ &\small =1+\frac1{\sqrt{\pi}} \int_{\frac1{2x}}^{0} s^{-\frac12}e^{-s}\,ds\\ &\small =1-\frac1{\sqrt{\pi}} \int_{0}^{\frac1{2x}} s^{-\frac12}e^{-s}\,ds \end{align*}\]  さらに \(\small s^{\frac12}=u\) とおくと, \(\small s=u^2,~ds=2udu\) であることから, \[\begin{align*} \small {\rm cdf\_}&\small {\rm levi}(x)\\ &\small =1-\frac2{\sqrt{\pi}} \int_0^{\frac1{\sqrt{2x}}} e^{-u^2}\,du\quad (1) \end{align*}\] となります.第2項は誤差関数です. 誤差関数は \[\small {\rm erf}(x)=\frac2{\sqrt{\pi}}\int_{0}^{x}e^{-t^2}\,dt\] で定義されるので,累積確率分布関数は \[\begin{align*} \small {\rm cdf\_}&\small {\rm levi}(x)\\ &\small =1-{\rm erf}\left(\frac1{\sqrt{2x}}\right) \end{align*}\] と表されます.
 乱数を発生させるには,この関数の逆関数を求めなければなりません. 上式から \[\small 1-{\rm erf}\left(\frac1{\sqrt{2y}}\right)=x\] とすると, \[\begin{gather*} \small {\rm erf}\left(\frac1{\sqrt{2y}}\right)=1-x\\ \small \frac1{\sqrt{2y}}={\rm erf}^{-1}(1-x)\\ \small \therefore\quad y=\frac1{2\left\{{\rm erf}^{-1}(1-x)\right\}^2} \end{gather*}\]  したがって,レビ分布の累積確率分布関数の逆関数は 次の式で表されます. \[\begin{align*} \small {\rm quant\_}&{\rm levi}(x)\\ &\small =\frac1{2\left\{{\rm erf}^{-1}(1-x)\right\}^2} \end{align*}\]  なお,「quantile_levi」とすべきところを簡略化しています. 誤差関数に対して \(\small {\rm erfc}(x)=1-{\rm erf}(x)\) を 相補誤差関数といいます. この関数を利用すると, \[\small {\rm quant\_levi}(x)=\frac1{2\left\{{\rm erfc}^{-1}(x)\right\}^2}\] と表すこともできます. 以下では \(\small {\rm erf}^{-1}(1-x)\) を利用した式で考えます. したがって,レビ分布にしたがう乱数は, \[\begin{align*} &\small {\rm random}\_{\rm levi}(1,0)\\ &\small ={\rm quant}\_{\rm levi}({\rm random}(1.0)) \end{align*}\] により発生させることができます.

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

グラフの確認
 次に,以上で定義した関数のグラフを確認してみましょう.

 (%i5)〜(%i7)では,それぞれ「pdf_levi」「cdf_levi」「quant_levi」 を定義しています.(%i8)(%i9)では, cdf_levi(x)とquant_levi(x)が互いに逆関数であることを 値の上から確認しています.
 グラフは,下図のようになります.

 ここで,確率密度関数に含まれる個々のグラフをみると, 下図のようになります.

 下図は quant_levi(x) の定義域を実数全体に拡張して cdf_levi(x) と同時表示し,直線 \(\small y=x\) に関して 対称であることを確認したものです.

 下図は,quant_levi(x) だけのグラフです.縦軸を対数軸にしました.

 具体的な値を求めると,「quant_levi(0.9)=63.33」 「quant_levi(0.99)=6365.86」 「quant_levi(0.999)=636619.4」などとなり, 極めて裾の重い分布であることが分かります.

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

乱数の分布
 以上の準備のもとで,レビ分布にしたがう標本の標本平均の分布の様子を 調べてみましょう.

 (%i13)は,累積確率分布関数の逆関数に一様乱数を利用して, レビ分布にしたがう乱数 \(\small n\) 個からなる リストを作成するコマンドです(参照). (%i14)は,\(\small n\) 個の平均からなる 要素数 \(\small N\) のリストを作成するコマンドです. 乱数がレビ分布にしたがっていることを確認するため, (%i16)では1個の標本平均からなる10000要素のリストを作成していますが, エラーが生じます. おそらくは,「quant_levi(x)」を定めている「inverse_erf(x)」, つまり誤差関数の逆関数の計算の際のエラーと思われ, 定義域での端点に近い箇所の値がうまく計算できなかった のだろうと思われます. 何度か試すとエラーが生じるような値が回避されます. この場合は(%i19)のように3回目でlevi1が生成されました.

 (%i20)(%i21)では,最小値と最大値を確認しています. レビ分布は裾が長いことから, (%i22)では区間 \(\small [0,100]\) を 200 等分したときの levi1 の度数分布を作成しています. (%i23)で \(\small x\) 軸の分点リストを定めます. (%i25)では,離散グラフを描画するときの書式をあらかじめ定めています. 度数を「階級幅×総度数」で割ってヒストグラムにしています.
 以上をもとにすると下図のようなグラフが得られます. 刻み幅を 0.5 にしているので最初の部分はちょっと切れますが, レビ分布の確率密度曲線と重なっていることが分かります.

 以下は,両対数グラフです.区間は [0,100] としました.


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

平均と分散
 レビ分布は,平均も分散も存在しません. 実際の計算で確かめてみましょう. \[\begin{align*} \small E(X) &\small =\int_{0}^{\infty}x\cdot{\rm pdf}\_{\rm levi}(x)\,dx\\ &\small =\int_{0}^{\infty}x\cdot\frac1{\sqrt{2\pi}} \frac1{x^{\frac32}}e^{-\frac1{2x}}\,dx\\ &\small =\int_{0}^{\infty}\frac1{\sqrt{2\pi}}x^{-\frac12}e^{-\frac1{2x}}\,dx\\ &\small \gt \int_{1}^{\infty}\frac1{\sqrt{2\pi}}x^{-\frac12}e^{-\frac12}\,dx\\ &\small =\frac1{\sqrt{2\pi e}}\int_{1}^{\infty}x^{-\frac12}\,dx\\ &\small =\frac1{\sqrt{2\pi e}}\big[2\sqrt{x}\big]_1^{\infty}=\infty \end{align*}\] ここでは,\(\small [1,\infty)\) では \(\small 1/2x\lt 1/2\) より \(\small -1/2x\gt -1/2\) であることを利用しています. したがって,平均が存在しないので,分散も存在しません.
 次に,実際の標本平均を調べて,「平均も分散も存在しない」ということを 確かめてみましょう.

 ここでは,10個の乱数の標本平均値を1000個持つリストを「levi10mean」, 10個の乱数の標本標準偏差を1000個持つリストを「levi10std」とします. 横軸を「回目」にして値をプロットすると下図のようになります. 値が大きく振れるので,いずれも縦軸を対数軸にしています. 上段が標本平均,下段が標本標準偏差ですが,いずれも値が大きく振れているのが 分かります. これは,標本数を増やしても同様です.つまり, 標本数を幾ら増やしても,標本平均や標本標準偏差を何回調べても, どのような値になりそうなのか全く分からないということです.

 具体的に,どのような値が乱数で生成されているのでしょうか.

 標本平均と標本標準偏差が大きな値になるリストを探すと,たとえば 標本平均が 994.1,標本標準偏差が \(\small 2.063\cdot 10^3\) になるものがあり, 次のような値でした.
\(\small 6.698\cdot 10^3, 5.162, 1.333, 486.5, 0.2898,\) \(\small 0.2355, 24.0, 2.712\cdot 10^3, 10.82, 2.99\)

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

平均の分布
 標本平均の値が大きく振れていることは分かりましたが, その値はどのような分布をしているのか2個平均と5個平均の場合で調べてみましょう. 他の分布の場合と同様に平均を要素とするリストを作成するには, 5個平均の場合は5万回の乱数を発生させることになります. 1万回でもエラーが生じます. 5回連続でエラーが出ないことは考えにくいので, 最初に,エラーの生じない1万要素のリストを生成しておくことにします. 1万個であれば,数回試せば生成されます.

 (%i31)〜(%i36)で1万個ずつ4万個の乱数を生成したことになります. 実際には(%i28)から生成を開始しているので, その間の番号が飛んでいる箇所はエラーが発生した箇所です. (%i37)(%i38)では, 作成したリストをまとめて2万個と5万個のリスト Llevi2, Llevi5 を作成しています. 「append」は,複数のリストを一つのリストにまとめるコマンドです. (%i39)(%i40)では,念のため要素数を確認しています.

 (%i41)では,2万個の乱数からなるリスト「Llevi2」の 最初から2個ずつの標本平均を計算し,それらを要素に持つリストを作成しています. (%i42)では,「Llevi5」をもとに5個ずつの標本平均を計算しています. いずれも,要素数は1万個のリストになります. 標本平均のリストをこのような形で作成したので, 実際には(%i14)のコマンドは不用になります.
 (%i43)(%i44)では,[0,100] を 200 等分した度数分布を作成し, (%i56)(%i57)でグラフを描くときの書式を定めています. 番号が飛んでいるのは,\(\small y\) 軸に指定するリスト名を 間違えていたことに気づいて修正したためです. \(\small x\) 軸の分点リストは,すでに(%i23)で作成済みです.
 以上をもとにグラフを描くと,次のようになります.

 以下は,両対数グラフです.範囲は [0,100] としています.

 いずれも,レビ分布の確率密度関数のグラフと類似したものになっています.

 確率分布の範疇に「安定分布」 と呼ばれるものがあり,レビ分布は安定分布です. そして,安定分布の性質[15]によれば, レビ分布の特性指数は \(\small \alpha=1/2\) であるので \(\small n^{1/\alpha}=n^2\) となり,

\[\small X_1+X_2+\cdots+X_n\overset{d}{=}n^2X\] となるはずなので,標本平均の分布については \[\small \frac1{n}(X_1+X_2+\cdots+X_n)\overset{d}{=}nX\] になるはずです. そこで,5個平均の分布について,このことを確かめてみましょう. 下記は,後でやり直しているので行番号がずれています.

(%i51)では,最初に生成したレビ分布にしたがう1個だけの平均からなる 1万個のリストを5倍しています.これを \(\small 5X\) とします. そして,同じような度数分布表を作って levidosu5b とし, 5個平均からなる1万個のリストの度数分布と重ねています. ほぼ同じ分布になることが分かります. コーシー分布の場合は, 何個の平均で考えても, その分布はもとの分布と一致しましたが,レビ分布の場合は \(\small n\) 個の平均の分布は,もとの分布を \(\small n\) 倍したものと 一致することが分かります. このような確率分布が,「安定分布」と呼ばれる確率分布です.

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


copyright