Maximaを活用した数学学習

in [数ナビの部屋]

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

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

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


[御案内] 「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 X\) の 確率密度関数を \(\small f(x)\) とすると, \(\small f(x)\) は次の性質を持ちます. 離散的な確率変数では,\(\small \int\) を \(\small \Sigma\) で 考えることになります. \[\small \begin{align*} &f(x)\geq 0\\ &\int_{-\infty}^{\infty}f(x)\,dx=1\\ &{\rm P}(a\leq X\leq b)=\int_a^bf(x)\,dx \end{align*}\]  確率密度関数 \(\small f(x)\) に対して, 累積分布関数 \(\small F(x)\) は次の式で定義されます. \[\small F(x)=\int_{-\infty}^{x}f(t)\,dt\]  微積分の基本定理により,両辺を \(\small x\) で微分すると, \[\small F'(x)=\frac{d}{dx}\int_{-\infty}^{x}f(t)\,dt=f(x)\] となります.
 統計では,累積分布関数の値を \(\small p\) とするとき, \[\small F(x)=\int_{-\infty}^xf(t)\,dt=p\] を満たす \(\small x\) を求める場面が頻出します. これは,累積分布関数の逆関数を考えることになります. このようなときのために, p分位数を求める関数が用意されています. たとえば,第一四分位数は, \[\small F(x)=\int_{-\infty}^{x}f(x)\,dx=1/4\] を満たす\(\small x\) のことです. 四分位数に限らず, 任意の \(\small p~(0\leq p\leq 1)\) に対して, \(\small p\) を与えて上式を満たす \(\small x\) の値を 求めることができます.
 確率変数 \(\small X\) の確率密度関数を \(\small f(x)\) とすると, 平均や分散は次の式で定義されます. \[\small \begin{align*} (平均)\quad &E[X]=\int_{-\infty}^{\infty}xf(x)\,dx\\ (分散)\quad &V[X]=\int_{-\infty}^{\infty}(x-E[X])^2f(x)\,dx \end{align*}\]  標準偏差は \(\small \sigma[X]=\sqrt{V[X]}\) です. 平均を \(\small \mu\) とするとき, 分散を \(\small (X-\mu)^2\) の平均とみて, \(\small V[X]=E[(X-\mu)^2]\) と表す場合もあります.
 多くの場合は以上のものだけで議論が進んで行きますが, 他に「歪度」や「尖度」というものもあります. それぞれ次の式で定義されます. \[\small \begin{align*} (歪度)\quad &E\bigg[\left(\frac{X-\mu}{\sigma}\right)^3\bigg]\\ (尖度)\quad &E\bigg[\left(\frac{X-\mu}{\sigma}\right)^4\bigg]-3 \end{align*} \]  歪度は分布の対称度を表す指標とされ, 歪度が \(\small 0\) のときは左右対称であり, 正のときは右側に裾が長く(左側にこぶがある), 負のときは左側に裾が長い(右側にこぶがある)分布になります. 尖度は分布の尖り具合を表す指標とされ, 尖度が \(\small 0\) のときは正規分布と同じような尖り具合で, 正のときは正規分布よりも尖っており, 負のときは正規分布よりも鈍い分布になります.

 Maximaでの記号は, 基本的に「(統計量)_(確率分布の名称)」で表され, 正規分布(normal distribution) N(\(\small m,s^2\)) の場合は次のように表されます. \(\small m\) は平均,\(\small s^2\) は分散です. したがって,\(\small s\) は標準偏差です.

確率密度関数:pdf_normal(x,m,s)
累積分布関数:cdf_normal(x,m,s)
p分位数:quantile_normal(p,m,s)
平均:mean_normal(m,s)
分散:var_normal(m,s)
標準偏差:std_normal(m,s)
歪度係数:skewness_normal(m,s)
尖度係数:kurtosis_normal(m,s)
乱数:random_normal(m,s)
 主な確率分布の名称は,次のように表されます.
二項分布:binomial
ポアソン分布:poisson
正規分布:normal
\(\small t\) 分布:student_t
\(\small \chi^2\) 分布:chi2
F分布:f
 繰り返しになりますが, Maximaは全部で20以上の確率分布について, 上記にあげたような確率密度関数や累積分布関数などが登録されているので, 通常現れるような確率分布には全て対応できると思われます. 詳細は,マニュアルの 「 52. distrib」を参照してください.


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

標準正規分布

 前述のMaximaのコマンドを利用すれば, 具体的な関数の式を知らなくても計算を進めていくことができますが, 以下では どのような関数で計算されているのかを確認しておきましょう. 確率分布への理解を深めるには,計算結果だけはなく, その値がどのようにして得られているのかを理解しておくことも必要です. 以下では標準正規分布 N(0,1) の場合を考えます( 参照).


関連する諸関数
 標準正規分布は,正規分布 N(m,s) で \(\small m=0, s=1\) の場合なので, \(\small f(x)={\rm pdf\_normal}(x,0,1)\), \(\small F(x)={\rm cdf\_normal}(x,0,1)\) です. そして,\(\small F(x)=p\) を満たす \(\small x\) は, \(\small x={\rm quantile\_normal}(p,0,1)\) と表されます. 具体的な関数は, \[\small f(x)=\frac1{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\] が確率密度関数であり, 累積分布関数はこの関数を積分することで得られます. \[\small F(x)=\frac1{\sqrt{2\pi}}\int_{-\infty}^{x}e^{-\frac{t^2}{2}}\,dt\]
 下図は,2つのパッケージ「distrib」と「descriptive」を読み込んでから 実行しています.

 (%i2)では \(\small -\infty\) からの広義積分を計算するため, 「integrate」ではなく「ldefint」で定義しています.「ldefint」は, 積分の両端の値を必要であれば極限値で計算します. また,(%i5)(%i6)のように \(\small x\) を含んだまま実行すると, 確率密度関数や分布関数の定義式が表示されます.(%o5)により, 確率密度関数が上述の \(\small f(x)\) と同じであることが分かります. (%o6)では「erf(x)」という関数を用いて表示されています. erf(x)は誤差関数と呼ばれ次の式で定義されます. \[\small {\rm erf}\,(x)=\frac{2}{\sqrt{\pi}}\int_{0}^xe^{-t^2}\,dt\]  この積分は初等積分可能ではないので, 既知の関数を利用して表すことはできません. 積分の式のまま定義するしかないのですが, 統計や実験系では重要な関数です. 下図は,誤差関数のグラフです.

 この誤差関数を利用すると分布関数 \(\small F(x)\) は次のように表され, (%o6) と同じものであることが分かります.下記の計算では, 途中で \(\small t/\sqrt{2}=s\) として置換積分を行っています. \[\small \begin{align*} &\small F(x)\\ &\small =\frac1{\sqrt{2\pi}}\int_{-\infty}^{x}e^{-\frac{t^2}{2}}\,dt\\ &\small =\frac1{\sqrt{2\pi}}\int_{-\infty}^{0}e^{-\frac{t^2}{2}}\,dt\\ &\small \qquad +\frac1{\sqrt{2\pi}}\int_{0}^{x}e^{-\frac{t^2}{2}}\,dt\\ &\small =\frac12+\frac1{2\sqrt{2}}\cdot\frac2{\sqrt{\pi}}\int_{0}^{x}e^{-\frac{t^2}{2}}\,dt\\ &\small =\frac12+\frac1{2\sqrt{2}}\cdot\frac2{\sqrt{\pi}} \int_0^{\frac{x}{\sqrt{2}}}e^{-s^2}\,\sqrt{2}\,ds\\ &\small =\frac12+\frac12{\rm erf} \left(\frac{x}{\sqrt{2}}\right) \end{align*}\]  これは,(%o6)で表示されている式と同じものです. 次に,p分位数を求める関数 (累積分布関数の逆関数) を表示させてみましょう.

 上図の(%o7)では関数の式がそのまま返されています. そこで,\(\small 0\lt x\lt 1\) であるということを 「assume」を利用して仮定すると,(%o9)にように式の形が表示されます. 「inverse_erf」は,誤差関数の逆関数ということです. 実際, \(\small F(x)=\frac12+\frac12{\rm erf} \left(\frac{x}{\sqrt{2}}\right)=p\) の逆関数は,次のように求められます.
 なお,最新版のMaximaでは,(%i8)の仮定をしなくても(%o9)が表示されます. \begin{gather*}\small \frac12+\frac12{\rm erf} \left(\frac{x}{\sqrt{2}}\right)=p\\ \small {\rm erf} \left(\frac{x}{\sqrt{2}}\right)=2p-1\\ \small \frac{x}{\sqrt{2}}={\rm erf}^{-1}(2p-1)\\ \small \therefore\quad x=\sqrt{2}\,{\rm erf}^{-1}(2p-1) \end{gather*}  「assume」を利用すると,いろいろな仮定を置くことができます. どのような仮定がなされているかは「facts()」により確認できます. 仮定を削除するには「forget」を利用します. (%i13)では,「facts()」により何も表示されないので, 仮定が削除されたことが分かります.


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

グラフの確認
 次に,グラフを確認しておきましょう.

 上図は,標準正規分布の確率密度曲線と累積分布曲線です.
 次に,累積分布曲線とp分位数を表す関数のグラフを描画してみましょう. 下図のように,直線 \(\small y=x\) に関して対称になり, 逆関数であることが分かります.

 ただし,quantile_normal(x,0,1) の定義域は \(\small 0\leq x\leq 1\) であり,定義域が全範囲である累積分布曲線と同列に扱うことができません. そこで,定義域の範囲外では 0 を取ることにして定義域を全範囲に拡張しておきます. それは,(%i19)のように if-then-else を利用すると定義できます. qn(x)は,次のような関数です. (%i20)の「same_xy」は,2つの座標軸のスケールを揃えるための指定です. デフォルトでは,ちょっと横長の画面に表示されます. \[\small \begin{cases} \quad 0 \hfill (x\lt 0)\\ {\rm quntile\_normal} (x,0,1)&\\ \hfill (0\leq x\leq 1)\\ \quad 0 \hfill (1\lt x)\end{cases}\]

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

幾つかの値の確認
 今度は,具体的な値を確認してみましょう. 統計解析では \(\small {\rm P} (|X|\leq \alpha)=0.95\) を満たす \(\small \alpha\) の値が頻出します. N(0,1)の確率密度曲線が左右対称であることから, このような \(\small \alpha\) は \(\small {\rm P} (X\leq \alpha)=0.975\) を満たすので, 要するに \(\small F(\alpha)=0.975\) となる値を求めれば良いわけです. それは,累積分布関数の逆関数なので,p分位数を表す関数を利用すると \(\small \alpha={\rm quantile\_normal} (0.975,0,1)\) により求められます.

 (%i21)の結果は(%o21)のように \(\small \sqrt{2}\) が残るので, (%i22)では直前の結果(%)を小数に直すために「numer」を付加することにより \(\small \alpha=1.96\) が得られています. \(\small \sqrt{2}\) が残るのは,この値が(%o9) により計算されているためと思われます.
 下図では,確率密度関数 \(\small f(x)\) の \(\small \alpha\) までの積分を実際に計算して0.975になることを 確認しています.「%」は,直前の結果である(%o22)のことであり, 「minf」は \(\small -\infty\) のことです.

 (%o23)では誤差関数を用いた結果が返されますが, 「numer」を付記して小数に直すと「0.975」が得られます. (%i23)の積分は,定積分を拡張した広義積分です. 「integrate」は広義積分を数値的に求めようとしますが, Maximaには広義積分を求めるコマンド「ldefint」も用意されています. 「ldefint」は, 積分の両端での値を必要であれば極限値の計算をして求めようとします. 下図のように, ldefintを利用することで広義積分の値が(%o23)よりもすっきりした形で 求められます.

 多くの統計の教科書では,標準正規分布N(0,1)に対して 「1.96」の値や \(\small {\rm P} (a\leq X\leq b)\) の値は正規分布表を 用いて計算するしかなく, その値がどのようにして計算されかは度外視しています. しかし,数式処理ソフトを利用すると, 結果を鵜呑みにすることなく, 定義式をもとに自分で直接計算して求めることができます. 以上の計算をしながら統計の学習を進めていけば, より理解が深まると思われます.

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

平均と分散
 次に,平均や分散の値を求めてみましょう.

 (%i29)(%i30)では,Maximaのコマンドを使って平均と分散を求めています. (%i31)(%i32)では,平均と分散の定義式により計算で求めています. 広義積分による計算なので「ldefint」を利用しています. (%o32)より,\(\small V[X]=E[X^2]-(E[X])^2=1\) が分かります. しかし,この値をみても「分かった!」感がしないかもしれません. そこで,定義式である積分の式を直接計算してみましょう.

 (%i33)は,N(0,1)の平均の定義式が \[\small \begin{align*} E[X]&=\int_{-\infty}^{\infty}tf(t)\,dt\\ &=\lim_{x\to\infty}\int_{-\infty}^{x}tf(t)\,dt \end{align*}\] と表されることから,\(\small (-\infty,x]\) の積分を求めて, その結果に対して極限値(limit)を求めています. 自分で計算するには,\(\small (-\infty,0]\) と \([0,\infty)\) に分けて 計算します. \[\begin{align*} \small \int_{-\infty}^{\infty}tf(t)\,dt &\small =\frac1{\sqrt{2\pi}}\int_{-\infty}^{0}te^{-\frac{t^2}{2}}\,dt\\ &\small \qquad +\frac1{\sqrt{2\pi}}\int_{0}^{\infty}te^{-\frac{t^2}{2}}\,dt\\ &\small =\frac1{\sqrt{2\pi}}\bigg[-e^{-\frac{t^2}{2}}\bigg]_{-\infty}^{0}\\ &\small \qquad +\frac1{\sqrt{2\pi}}\bigg[-e^{-\frac{t^2}{2}}\bigg]_{0}^{\infty}\\ &\small =-\frac1{\sqrt{2\pi}}+\frac1{\sqrt{2\pi}}\\ &\small =0 \end{align*}\] となるので 0 になります.被積分関数が奇関数でも, それぞれの広義積分が存在しない場合もあるので注意が必要です.
 次に,分散について同様のことを行ってみましょう.

 (%i35)では,平均E[X]のときと同様のことを行っています. (%i37)では,不定積分を求めて, (%i38)(%i39)で個別に極限値の値を求めています. \(\small \frac12-(-\frac12)=1\) となります. \(\small e^{-\frac{x^2}{2}}\) は初等積分可能ではありませんが, \(\small xe^{-\frac{x^2}{2}}, x^2e^{-\frac{x^2}{2}}\) は 不定積分が存在します.ついでに, \(\small x^2f(x)\) のグラフも確認しておくとよいでしょう. 下図のグラフになります.


 以上のように,数式処理ソフト「Maxima」を利用すると, 統計に関する数表が手元に無くても統計分析を行うことができます. それはエクセルやRでも同様ですが, 「数式処理」機能を利用することで, いろいろな統計量の定義に基づいた積分計算を行うことができます. その計算も交えながら学習していくと, 統計に関する理解がより深まるのではないかと思われます.

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

平均の分布
 標準正規分布 N(0,1) から無作為抽出した大きさ \(\small n\) の標本の 平均は,平均が 0, 分散は \(\small 1/n\) の正規分布にしたがいます. 乱数を利用して,このことを実際に確かめてみましょう.

 (%i9)は,N(0,1)にしたがう乱数 \(\small n\) 個からなるリストを 作成するコマンドです. (%i10)では norlist(n) の平均からなる大きさ \(\small N\) のリストを, (%i11)では 標準偏差からなる大きさ \(\small N\) のリストを作成する コマンドです. (%i15)(%i16)では,5個の平均と標準偏差からなる大きさ1000のリストを 作成しています. 以上をもとに,「mlist」と「slist」の値の変化をみると 次のようになります.

 上図は,平均を計算するごとにどのような値になったかを グラフ化したものです.横軸は何回目に計算した平均であるかを表します. 0 を中心に多くは \(\small -1\sim 1\) の間の値になっているようです.
 下図は,標準偏差についてみたものです. 標準偏差は,1 を中心に多くは \(\small 0.5\sim 1.5\) の間の 値になっているようです.

 次に,5個平均からなるリスト「mlist」の度数分布を求めて ヒストグラムを描画してみましょう.

 (%i22)では,区間 \(\small [-5,5]\) を50等分したときの度数を得ています. (%i23)は,対応する横軸の分点リストです. (%i24)では,グラフ描画のための諸式です. 度数をヒストグラムの値に直しています.
 以上のもとでグラフを描画したのが下図になります. 5個平均の分布はN(0,1/5)にしたがうので, その確率密度関数のグラフも描画しました. 重なっているのが分かります.この場合の確率密度関数は, \(\small {\rm pdf}\_{\normal}(x,0,1/\sqrt{5})\) です.


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

無限分割可能分布
 \(\small X_1,X_2,\ldots,X_n\) をN(0,1)にしたがう互いに独立な 確率変数とすると,\(\small X_1+X_2+\cdots+X_n\) は 正規分布N(0,1/n)にしたがいます.これを逆に考えると, \(\small X_1,X_2,\ldots,X_n\) をN(0,1/n)にしたがう互いに独立な 確率変数とすると,\(\small X_1+X_2+\cdots+X_n\) は 標準正規分布N(0,1)にしたがうことになります. つまり,標準正規分布にしたがう確率変数を \(\small Z\) とすると, \(\small Z\) は,正規分布N(0,1/n)にしたがう確率変数の \(\small n\) 個の 和で表されることになります.\(\small n\) は任意の自然数です.
 このように,一つの確率変数 \(\small Z\) を, 別な確率変数の和で表すことができるとき,確率変数 \(\small Z\) の したがう確率分布を「無限分割可能分布」といいます.
 このことを実際に確かめてみましょう.

 (%i26)で N(\(\small 0, 1/\sqrt{5}\)) にしたがう乱数のリストを作り, (%i27)ではその和からなるリストを作ります. (%i28)で,5個の和からなる大きさ1000のリストを「nor5sum」として, (%i29)で「nor5sum」の度数分布を得ています.
 グラフを描くと下図のようになり,標準正規分布のグラフと 重なっているのが分かります.


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

他の確率分布との関係
 いろいろな確率分布の中で, 標準正規分布N(0,1)は基本的で重要な確率分布です. それは,他のいろいろな確率分布と密接に関わっているからです.
 以下では,\(\small Z,W\) や \(\small Z_1,Z_2,\ldots,Z_n\) は標準正規分布にしたがうものとします.
  1. \(\small \mu\) を実数,\(\small \sigma\) を正の実数とすると, \(\small X=\mu+\sigma Z\) は正規分布 N(\(\small \mu,\sigma^2\)) に したがう.
  2. \(\small X=Z^2\) とおくと,\(\small X\) は自由度1の カイ二乗分布にしたがう.
  3. \(\small X=Z_1^2+Z_2^2+\cdots+Z_n^2\) とおくと,\(\small X\) は自由度 \(\small n\) のカイ二乗分布にしたがう.
  4. \(\small X=1/Z^2\) とおくと,\(\small X\) は レビ分布にしたがう.
  5. \(\small X=Z/W\) とおくと,\(\small X\) は コーシー分布にしたがう.
  6. \(\small X=e^Z\) とおくと,\(\small X\) は 対数正規分布にしたがう.
  7. \(\small X=|Z|\) とおくと,\(\small X\) は Folded Normal Distributionにしたがう.
  8. \(\small X=\sqrt{Z_1^2+Z_2^2}\) とおくと,\(\small X\) は レイリー分布にしたがう.
  9. \(\small X=\sqrt{Z_1^2+Z_2^2+Z_3^3}\) とおくと,\(\small X\) は マックスウェル分布にしたがう.

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

中心極限定理
 中心極限定理は,次のような定理です.
 確率変数 \(\small X_1,X_2,\ldots,X_n\) が 平均 \(\small \mu\),分散 \(\small \sigma^2\) を持つ 同一の確率分布にしたがうものとする.このとき, その平均を \(\small \bar{X}=\frac1{n}(X_1+X_2+\cdots+X_n)\) とすると, 次のことが成り立つ. \[\small \begin{align*} &\lim_{n\to\infty} {\rm P} \left(\frac{\bar{X}-\mu}{\sigma/\sqrt{n}}\leq a\right)\\ &=\frac1{\sqrt{2\pi}}\int_{-\infty}^{a} e^{-\frac{x^2}{2}}\,dx\end{align*}\]  これは,要するに, 「平均や分散が存在するような分布であれば, それがどのような分布であっても, 標本数を大きくしていけば, 標本平均を標準化した値の分布は標準正規分布N(0,1)に近づいていく」 という内容の定理です. 統計での多くの計算が結局は正規分布に帰着するのは, この定理によるところが大きいといえるでしょう. 統計を理解する上では, この定理の意味を正しく把握することが重要です. そのためには「どのような分布」でも, 標本数を大きくすると「平均の分布が正規分布に近づいていく」 ことを実感として理解することが必要です.
 ここでは,一様分布について以下のことを行ってみましょう.
  1. 個々の確率分布について,その分布にしたがう乱数を発生させる.
  2. 発生させた乱数の平均を求める.
  3. 以上のことを多数回行って,平均からなるリストを作成する.
  4. 作成されたリストから度数分布を求める.
  5. 平均を取る標本数を増やすと標準正規分布に近づくことを確認する.
 一様分布にしたがう乱数は「random」により発生させることができるので, 以下のような手順で確かめることができます.
  1. 区間 \(\small [0,1)\) 内の乱数は「random(1.0)」により得られるので, 乱数の \(\small n\) 個のリストを 「makelist(random(1.0),i,1,n)」により作成する. それを,「rlist(n)」とする.
  2. 作成したリスト「rlist(n)」の平均を「mean」により 求めて「heikin(n)」とする.
  3. \(\small n\) 個の平均からなる大きさ \(\small N\) のリストを 「makelist」で作成して「hlist(n,N)」とする.そして, そのリストを標準化して「shlist(n,N)」とする.
  4. 標準化されたリスト「shlist(n,N)」の度数分布を 「continuous_freq」で求める. 表示される結果の第2成分が度数分布である.
  5. 区間 \(\small [-5,5]\) を20等分したときの ヒストグラムを作成して,標準正規分布のグラフと比較する.
 一般に,区間 \(\small [a,b]\) における一様分布の 確率密度関数 \(\small f(x)\) は \[\small \begin{cases} ~0 & (x\lt a)\\ \frac1{b-a}& (a\leq x\leq b)\\ ~0& (b\lt x) \end{cases}\] で表されます.確率密度関数 \(\small f(x)\) は, 単位ステップ関数を利用すると \[\small \frac{{\rm unit\_step}(x-a)-{\rm unit\_step}(x-b)}{b-a}\] により表されます.区間 \(\small [0,1]\) では, \(\small f(x)={\rm unit\_step}(x)-{\rm unit\_step}(x-1)\) です. この関数の積分を考える必要があるので, 「abs_integrate」を読み込んでおきます (参照).

 上図では,入出力番号をリセットして最初からやり直しています. 「distrib」と「descriptive」を読み込んでから (%i5)で \(\small f(x)\) を定義して, この関数の積分を考えるために(%i6)で「abs_integrate」を読み込みます. (%i7)で \(\small f(x)\) が確率密度関数の条件を満たすことを 確認して, (%i8)では平均 \(\small \int_{-\infty}^{\infty}xf(x)\,dx\) の値を, (%i9)では \(\small \int_{-\infty}^{\infty}x^2f(x)\,dx\) の値を, それぞれ広義積分の値を求めることのできるコマンド「ldefint」で 計算します. (%i10)では,\(\small V[X]=E[X^2]-(E[X])^2\) という関係から 分散を求めています.

 上図では,(%i11)で \(\small n\) 個からなる乱数のリストを 作成するコマンド「rlist」を定義し, (%i12)でリストが作成されることを確認しています. (%i13)は,rlist(n) の平均で構成され, \(\small N\) の要素からなるリストを 作成するコマンド「hlist」を定義しています. (%i14)でリストが作成されることを確認しています. 個々の値は, それぞれ3個の乱数の平均を計算したもの(のはず)です. (%i28)は,平均からなるリスト「hlist(n)」を標準化しています. 「float」は,結果を浮動小数点で求めるためのものです. これをつけないと,\(\small \sqrt{12}=2\sqrt{3}\) であることから, 小数と \(\small \sqrt{3}\) が混在した形で結果が表示されます.

 (%i30)は,標準化されたリスト「shlist」の \(\small [-5,5]\) に おける度数分布を求めています.0.5 刻みにするため, 区間を20等分しています. 範囲はエラーが生じないよう大きめに取りました (参照). (%i31)では,\(\small N=100\) の場合を確認したものです. (%i32)では,\(\small x\) 軸の区分点のリストを作成しています. 0.5刻みの小区間の中点になるように設定しています. (%i33)は,plot2dでの離散データ部分の指定を, 事前に設定したものです.柱の面積が相対度数と一致するように \(\small 0.5*N\) で割っています.
 以上の準備のもとで,(%i38)では \(\small n=2, 10\) の場合の 平均の値を1000個用意した場合の分布図です.

 なお,(%i28)で標準化を行わないで度数分布を表示させると, \(\small n\) 個の平均の分布は, \(\small n\) の値が大きくなるにつれ元の平均 (1/2) の近くに集まってくる 様子を確認することができます.下図は,\(\small n=2, 5\) の場合に, 区間 \(\small [0,1]\) を20等分して相対度数分布を描画したものです. なお,(%o8)にあるように,この分布の平均は 1/2 です. Lx,dosu(n,N),Gr(n,N)を,それぞれ次のように修正しています. plot2dでの範囲は,[x,0,1], [y,0,0.2] としました.


 他の確率分布にしたがう場合も, 同じ手順で標準化した場合のグラフや, 標準化しない場合のグラフを確認することができます. 標準化の場合の修正箇所は下記の通りです.
(%i5)〜(%i10)
  • 確率密度関数 \(\small f(x)\) を修正する.
  • 確率密度関数は「pdf_(name)」の形で表される. 個々の確率分布の名称は 「 52. distrib」を参照.
  • 平均は「mean_(name)」,標準偏差は「std_(name)」により求められる. これらの値は,平均を標準化するときに必要となる.
(%i11)〜(%i29)
  • (%i11)の「random」を「random_(name)」に変更する. その引数は個々の確率分布に合わせて指定する.
  • (%i28)の標準化の式で平均と標準偏差の箇所を修正する.
(%i30)〜(%i38)
  • (%i38)の「Gr(n,N)」の \(\small n, N\) の値は, 自分で表示させたい値を指定する.

 中心極限定理については,下記も参照してください.

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

カイ二乗分布

 中心極限定理が成り立つような確率分布では, その分布から抽出した標本平均を標準化すると標準正規分布に近づいていきますが, 分散の分布はどのようになっているのでしょうか. この分散の分布を調べるときに必要になるのが \(\small \chi^2\) 分布です 参照).


正規乱数のニ乗和
 分散がどのような分布になるのか, 標準正規分布N(0,1)から \(\small n\) 個の標本を抽出して試してみましょう. \(\small X_1, X_2, \ldots, X_n\) を抽出したとすると, \(\small (平均) \mu=E[X]=0\) となるので, \(\small V[X]=\frac1{n}\sum_{k=1}^{n}X_i^2\) を調べることになります. まず,\(\small X_1^2+X_2^2+\cdots+X_n^2\) について調べてみましょう. この式がしたがう分布が自由度 \(\small n\) のカイニ乗分布です. 「desciptive」と「distrib」は読み込み済みとします.

 上図は,パッケージの「ditrib」「descriptive」を読み込み済みとします. (%i12)では,N(0,1)にしたがう正規乱数の2乗からなる, 大きさ \(\small n\) のリストを作成しています. (%i13)では,そのリストの要素の合計を「平均×個数」により計算して, その合計値からなる大きさ \(\small N\) 個のリストを作成しています. マニュアルを見てもリスト内の要素の合計を求めるコマンドが見当たらないので, それを「平均×個数」で求めました. (%i14)では,平方和の合計値からなるリストの要素の度数分布を, \(\small [0, 100]\) を200等分して求めています. 範囲が大きすぎるように思えますが, 正規乱数の箇所で述べたような理由で, エラーが生じないようにするための処置です. (%i15)では \(\small x\) 軸の区分点のリストを作成し, (%i16)ではplot2d内の離散データの書式を定めています.
 以上のもとで,正規乱数の2乗和の度数分布は次のように描画されます. \(\small n=1, 2, 4\) のときの1000個のデータの分布です.

 下図は,\(\small \chi^2\) 分布の確率密度関数(pdf_chi2(x,n))のグラフです.


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

確率密度関数
 ここで,\(\small \chi^2\) 分布の確率密度関数の定義式を確認してみましょう.

 (%o21)をみると,\(\small \chi^2\) 分布の確率密度関数は ガンマ分布の確率密度関数を用いて表されています. しかし,ガンマ分布の確率密度関数を見ようとしてもコマンドしか返されません. 周知のように, \(\small \chi^2\) 分布の自由度 \(\small n\) の確率密度関数は 次のようなガンマ関数を用いて表されます. \[\small \frac1{2^{\frac{n}{2}}\Gamma\left(\frac{n}{2}\right)} x^{\frac{n}{2}-1}e^{-\frac{x}{2}}\quad (x\gt 0)\]  通常は,この関数が天下りに提示されてそれを信じるしかないのですが, 最も簡単な \(\small n=1\) の場合であれば, 何故このような関数になるかをある程度理解することができます.

 確率変数 \(\small X\) が標準正規分布N(0,1)にしたがえば, その確率密度関数は次の式で表されます. \[\small f(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\] 自由度が \(\small n=1\) の場合は \(\small X^2\) の確率密度関数を 考えることになりますが, それは \(\small Y=X^2\) とおくとき,累積分布を考えると \[\small P(Y\leq x)=\int_{0}^{x}g(y)\,dy\] を満たすような関数 \(\small g(y)\) を求めることになります. \[\begin{align} \small P(Y\leq x) &\small =P(X^2\leq x)\\ &\small =P(-\sqrt{x}\leq X\leq \sqrt{x})\\ &\small =2P(0\leq X\leq \sqrt{x})\\ &\small =2P(X\leq \sqrt{x})-1 \end{align}\] これは,積分で表すと \[\small \int_{0}^{x}g(y)\,dy =2\int_{-\infty}^{\sqrt{x}}f(t)\,dt-1\] ということなので,両辺を \(\small x\) で微分すると 微積分の基本定理と合成関数の微分法により \[\begin{align*} \small g(x) &\small =2f(\sqrt{x})(\sqrt{x})'\\ &\small =2\cdot\frac1{\sqrt{2\pi}}e^{-\frac{x}{2}}\cdot\frac1{2\sqrt{x}}\\ &\small =\frac1{\sqrt{2\pi}}x^{-\frac12}e^{-\frac{x}{2}} \end{align*}\] となります. ここで,\(\small \sqrt{2}=2^{\frac12}\) であり, ガンマ関数の性質から \(\small \Gamma\left(\frac12\right)=\sqrt{\pi}\) であることを利用すると, この関数は次のように表すことができます. \[\small g(x)=\frac{1}{2^{\frac12}\Gamma\left(\frac12\right)}x^{\frac12-1}e^{-\frac{x}{2}}\]  以上のことを \(\small X_1^2+X_2^2+\cdots+X_n^2\) で考えて得られるのが, 前述した \(\small \chi^2\) 分布の確率密度関数です.

 周知のように,\(\small Z_1, Z_2, \ldots, Z_n\) がN(0,1)にしたがえば, その平方和 \(\small Z_1^2+Z_2^2+\cdots+Z_n^2\) は自由度 \(\small n\) の \(\small \chi^2\) 分布にしたがいます.このことを利用すると, \(\small X_1, X_2, \ldots, X_n\) がN(\(\mu\),\(\small \sigma^2\))にしたがえば, それを標準化した \[\begin{align*} \small \left(\frac{X_1-\mu}{\sigma}\right)^2 &\small +\cdots+\left(\frac{X_n-\mu}{\sigma}\right)^2\\ &\small =\frac1{\sigma^2}\sum_{k=1}^{n}(X_i-\mu)^2\end{align*}\] は自由度 \(\small n\) の \(\small \chi^2\) 分布にしたがいます. 母平均 \(\small \mu\) が不明のとき, それを標本平均 \(\small \bar{X}\) で置きかえた \[\small \frac1{\sigma^2}\sum_{k=1}^{n}(X_i-\bar{X})^2\] は自由度 \(\small n-1\) の \(\small \chi^2\) 分布にしたがうことが 知られています.


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

平均と分散
 前項で述べたように,自由度 \(\small n\) のカイ二乗分布の 確率密度関数 \(\small f(x)\) は, \[\small \frac1{2^{\frac{n}{2}}\Gamma\left(\frac{n}{2}\right)} x^{\frac{n}{2}-1}e^{-\frac{x}{2}}\quad (x\gt 0)\] です. (%o21)で示されるガンマ分布の確率密度関数は, Maximaのソースファイルが置かれているフォルダー内にある「distrib.mac」の中で, 次の式で記述されています. \[\small \frac{x^{a-1}e^{-x/b}}{b^a\Gamma(a)}\] したがって,カイ自乗分布の確率密度関数は,ガンマ分布の確率密度関数の \(\small a=\frac{n}{2}, b=2\) の場合に相当します. つまり,この式を \(\small h(x,a,b)\) で表すことにすると, \(\small g(x)=h(x,\frac{n}{2},2)\) ということです.
 そこで,最初にガンマ分布について考えてみましょう. 簡単のため \(\small b=1\) の場合を考えると, \[\small h(x,a,1)=\frac1{\Gamma(a)}x^{a-1}e^{-x}\] であり,ガンマ関数の定義式の被積分関数が含まれています. そのことを利用すると,\(\small b=1\) のときの平均は次のようになります. \[\begin{align*} \small E(X) &\small =\int_0^{\infty}xh(x,a,1)\,dx\\ &\small =\frac1{\Gamma(a)}\int_0^{\infty}x^{(a+1)-1}e^{-x}\,dx\\ &\small =\frac{\Gamma(a+1)}{\Gamma(a)}\\ &\small =\frac{a\Gamma(a)}{\Gamma(a)}=a \end{align*}\] ここでは, ガンマ関数 の \(\small \Gamma(x+1)=x\Gamma(x)\) という性質を利用しています. 同様にすると, \[\begin{align*} \small E(X^2) &\small =\int_0^{\infty}x^2h(x,a,1)\,dx\\ &\small =\frac1{\Gamma(a)}\int_0^{\infty}x^{(a+2)-1}e^{-x}\,dx\\ &\small =\frac{\Gamma(a+2)}{\Gamma(a)}\\ &\small =\frac{(a+1)a\Gamma(a)}{\Gamma(a)}\\ &\small =a(a+1) \end{align*}\] となるので,\(\small b=1\) のときの分散は次のようになります. \[\begin{align*} \small V(X) &\small =E(X^2)-E(x)^2\\ &\small =a(a+1)-a^2=a \end{align*}\]  次に,一般の \(\small h(x,a,b)\) の場合を考えると, \[\small h(x,a,b)=\frac1{b^a\Gamma(a)}x^{a-1}e^{-\frac{x}{b}}\] において,\(\small \frac{x}{b}=t\) とおくと \(\small x=bt\) である ことから, \[\begin{align} \small h(x,a,b) &\small =\frac1{b^a\Gamma(a)}(bt)^{a-1}e^{-t}\\ &\small =\frac1{b}\cdot\frac1{\Gamma(a)}t^{a-1}e^{-t}\\ &\small =\frac1{b}h(t,a,1) \end{align}\] となり,\(\small b=1\) の場合の式で表すことができます. したがって,この \(\small x=bt\) という置換を行うと, 一般の場合の平均は \[\begin{align*} \small E(X) &\small =\int_0^{\infty}xh(x,a,b)\,dx\\ &\small =\int_0^{\infty}bt\cdot \frac1{b}h(t,a,1)\,bdt\\ &\small =b\int_0^{\infty}t\cdot h(t,a,1)\,dt=ab \end{align*}\] 同様にすると \[\begin{align*} \small E(X^2) &\small =\int_0^{\infty}x^2h(x,a,b)\,dx\\ &\small =\int_0^{\infty}(bt)^2\cdot \frac1{b}h(t,a,1)\,bdt\\ &\small =b^2\int_0^{\infty}t^2\cdot h(t,a,1)\,dt\\ &\small =a(a+1)b^2 \end{align*}\] となるので,一般のガンマ分布の分散は \[\begin{align} \small V(X) &\small =E(X^2)-E(X)^2\\ &\small =a(a+1)b^2-(ab)^2\\ &\small =ab^2 \end{align}\] となります.ガンマ分布で \(\small a=\frac{n}{2}, b=2\) の場合が カイ二乗分布なので,自由度が \(\small n\) のカイ二乗分布の平均と分散は, 結局は次のような簡単な式で表されます. \[\small E(X)=n,\quad V(X)=2n\]

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

無限分割可能分布
 \(\small Z_i~(1\leq i\leq n)\) をN(0,1)にしたがう互いに独立な 確率変数とすると, \(\small Z_1^2+Z_2^2+\cdots+Z_n^2\) は自由度 \(\small n\) の カイ自乗分布にしたがいます.一般に, \(\small X_i~(1\leq i\leq n)\) が自由度 \(\small k\) のカイ二乗分布にしたがう互いに独立な確率変数とすると, \(\small X_1+X_2+\cdots+X_n\) は自由度 \(\small nk\) の カイニ乗分布にしたがいます. これを逆にいうと,\(\small X\) が自由度 \(\small N\) の カイニ乗分布にしたがい \(\small N=nk\) と因数分解できれば, \(\small X\) は自由度 \(\small k\) のカイ二乗分布にしたがう 互いに独立な確率変数 \(\small X_i~(1\leq i\leq n)\) に対して, \[\small X=X_1+X_2+\cdots+X_n\] と \(\small n\) 個の確率変数の和に分解できることになります.
 以上のことを,\(\small N=6=2\cdot 3\) の場合に 実際に確かめてみましょう. 前述のことから,\(\small X\) が自由度6のカイ二乗分布にしたがえば, \(\small X\) はそれぞれ自由度が2の3個の確率変数と, 自由度が3の2個の確率変数の和で表せることになります.

 (%i30)は,自由度が「fd」のカイ二乗分布にしたがう乱数「n」個の 和からなる,大きさ「N」のリストを作成するコマンドです. (%i31)では,そのコマンドを利用して,自由度が2のカイ二乗分布に したがう乱数3個の和からなる大きさ1万のリストを作成してます. (%i32)では,自由度が3のカイ二乗分布に したがう乱数2個の和からなる大きさ1万のリストを作成してます.

 (%i33)(%i34)では,作成した2つのリストから,区間 \(\small [0,100]\) を 500等分したときの度数分布を得ています. (%i36)(%i37)で,度数分布をヒストグラムの数値に直します.
 以上もとで,グラフを描画したのが下図です. 2つのヒストグラムが, 自由度6のカイ二乗分布の確率密度関数のグラフと 重なっているのが分かります.


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

t分布
 正規分布N(\(\small \mu,\sigma^2)\)から抽出した大きさ \(\small n\) の 標本を \(\small X_1, X_2, \ldots, X_n\) とすると, 周知のように標本平均 \(\small \bar{X}=\frac1{n}(X_1+X_2+\cdots+X_n)\) は 正規分布N(\(\small \mu,\sigma^2/n\))にしたがいます. 最初に,このことを確認してみましょう. 以下では,「distrib」と「descriptive」は読み込み済みとします 参照).

 (%i4)では \(\small \mu=50, \sigma=20\) の正規分布にしたがう乱数 からなる大きさ \(\small n\) のリストを作成しています. (%i5)は大きさ \(\small n\) のリストの平均からなる, 大きさ \(\small N\) のリストを作成しています. (%i11)では,区間 \(\small [0,100]\) を50等分したときの度数分布を 得ています. (%i12)は,小区間の中点からなるリストを作成しています. (%i13)は,plot2dで離散データのグラフを描画するときの書式が長いので, あらかじめ定めています.
 以上の準備のもとに,(%i20)で \(\small n=2, 10, 40\) の場合の 平均からなる500個の値を, 2点刻みの度数分布にしてグラフ化したものが下図です.

 平均を取る個数が増えるほど,平均の周りに近づいていきます. \(\small n\) 個の平均では分散が \(\small \sigma^2/n\) になるので, その理論曲線も描画してみましょう.

 (%i21)では,正規分布の確率密度関数を \(\small f(x,m,s)\) として, (%i22)で \(\small f(x,50,20)\) の式を表示させています. 正規分布の定義式から \(\small f(x,50,20)\) は次のようになり, (%o22)に表示されている式と一致します. \[\small \begin{align*} &\small f(x,50,20)\\ &\small =\frac1{\sqrt{2\pi}\cdot 20}e^{-\frac{(x-50)^2}{2\cdot 20^2}}\\ &\small =\frac1{5\cdot 4\sqrt{2}\sqrt{\pi}}e^{-\frac{(x-50)^2}{800}}\\ &\small =\frac1{5\cdot 2^{\frac52}\sqrt{\pi}}e^{-\frac{(x-50)^2}{800}} \end{align*}\]  大きさ \(\small n\) のときの標準偏差は \(\small 20/\sqrt{n}\) であるので, (%i23)により対応する3つの場合のグラフを描画すると, 下図が得られます.


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

 以上の分布は,母集団となる正規分布の分散 \(\small \sigma^2\) の 値が分かっている場合です.その場合は,標準化を行った \[\small Z=\frac{\bar{X}-\mu}{\sigma/\sqrt{n}}\] は,標準正規分布N(0,1)にしたがいます.  母分散 \(\small \sigma^2\) が不明のときは標本の値から推定するしかありません. それは,周知のように不偏分散 \(\small U^2\) で推定することができます. その場合に,標準化に相当する式は \[\small T=\frac{\bar{X}-\mu}{U/\sqrt{n}}\] となります.この値はどのような分布をするのでしょうか. \(\small T\) の分子・分母を \(\small \sigma/\sqrt{n}\) で割ると, \[\begin{align*}\small T&=\frac{\frac{\bar{X}-\mu}{\sigma/\sqrt{n}}} {\frac{U/\sqrt{n}}{\sigma/\sqrt{n}}}\\ &=\frac{Z}{U/\sigma} \end{align*}\] となります.さらに,分母は, \[\small\begin{align*} &\frac{U}{\sigma}=\sqrt{\frac{U^2}{\sigma^2}}\\ &=\sqrt{\frac1{n-1}\cdot\frac1{\sigma^2}\sum_{k=1}^{n}(X_i-\bar{X})^2} \end{align*}\] となります. ここで,\(\small X_1,X_2,\ldots,X_n\) は N(\(\small \mu,\sigma^2\))から抽出した標本とします. その場合,\(\small (X_i-\mu)/\sigma\) はN(0,1)にしたがうので, \(\small \frac1{\sigma^2}\sum_{k=1}^n(X_i-\mu)^2\) は 自由度 \(\small n\) の \(\small \chi^2\) 分布にしたがいますが, \(\small \mu\) を \(\small \bar{X}\) で置きかえた \[\small \frac1{\sigma^2}\sum_{k=1}^n(X_i-\bar{X})^2\] は,自由度 \(\small n-1\) の \(\small \chi^2\) 分布にしたがうことが 知られています. つまり,\(\small U/\sigma\) は, 自由度 \(\small n-1\) の \(\small \chi^2\) 分布にしたがう統計量を, その統計量の自由度である \(\small n-1\) で割って平方根を取った値です. したがって,\(\small T\) がどのような分布にしたがうかを調べるには, 標準正規分布N(0,1)にしたがう統計量を \(\small Z\), 自由度 \(\small n\) の \(\small \chi^2\) 分布にしたがう統計量を \(\small Y\) とするとき, \[\small T=\frac{Z}{\sqrt{\frac{Y}{n}}}\] がどのような分布にしたがうかを調べればよいことになります. この式が自由度 \(\small n\) の \(\small t\) 分布の定義式です.
 以上のことを実際に行ったのが下記です.

 (%i27)で \(\small T\) の値の大きさ \(\small n\) のリストを作ります. (%i45)では,\(\small [-20,20]\) を200等分した度数分布を求めます. 0.2刻みに区分したので,(%i46)は小区間の中点からなるリストです. (%i51)により,自由度が \(\small 5, 10, 20\) の場合に, 1万個の乱数をもとに0.2刻みにした場合のグラフが下記です.

 \(\small t\) 分布の確率密度関数も, 下記のようにガンマ関数を用いて表されます. Maximaでのコマンドは「pdf_student_t(x,n)」です. \[\small \frac{\Gamma\left(\frac{n+1}{2}\right)} {\sqrt{n\pi}\Gamma\left(\frac{n}{2}\right)} \left(1+\frac{x^2}{n}\right)^{-\frac{n+1}{2}} \quad (n\geq 1)\]  下図は,自由度が10の場合の1万個の分布と重ね合わせたものです. スペルが長いので,確率密度関数をf(x,n)として再定義しています. 描画するたびに新たな乱数が発生するので,上図の df=10 の 場合のグラフとは異なります.

 下図は,標準正規分布と比較したものです. 自由度が上がるごとにN(0,1)に近づいていきます. 自由度が30以上では,N(0,1)のグラフとほとんど一致します。 したがって,\(\small t\) 分布は標本数が30未満で少ない場合に必要な分布です.


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

コーシー分布

 中心極限定理は, 平均と分散が存在するような同一分布から標本を抽出した場合です. それらが存在しない分布とは,どのような分布なのでしょうか. それは,要するに,平均と分散の定義式 である広義積分が存在しない場合の分布です. そのような分布として,コーシー分布が知られています. 平均と分散のみならず,歪度も尖度も存在しません 参照).


確率密度関数
 コーシー分布の確率密度関数は,次の式で表されます. この式で表されるコーシー分布を Cauchy(a,b) で表すことにします. \[\small f(x,a,b)=\frac{1}{\pi}\cdot\frac{b}{(x-a)^2+b^2}\]

 (%i3)では,コーシー分布の確率密度関数を \(\small f(x,a,b)\) として 定義して,具体的な関数として \(\small f(x,2,3)\) を表示させています. \(\small a=0\) の場合で \(\small b\) の値を変えたグラフは, (%i7)により下図のようになります. \(\small b\) の値が大きくなるほどグラフは「のっぺり」としたものになり、 小さくなるほど尖度が大きくなります.

 このような分布は,どのような場面で生じてくるのでしょうか? コーシー分布は \(\small t\) 分布とも関わっており, Cauchy(0,1) は自由度 1 の \(\small t\) 分布と一致します. 実際,自由度 \(\small n\) の \(\small t\) 分布の 確率密度関数で \(\small n=1\) のときは \[\small \frac{\Gamma(1)} {\sqrt{\pi}\Gamma\left(\frac{1}{2}\right)} \left(1+x^2\right)^{-1}\] であり,ガンマ関数の性質から \(\small \Gamma(1)=1, \Gamma\left(\frac12\right)=\sqrt{\pi}\) であることに注意すると,この関数は Cauchy(0,1)の確率密度関数 \[\small \frac1{\pi}\cdot\frac1{1+x^2}\] と一致します。

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

N(0,1)の比
 Cauchy(0,1)は標準正規分布N(0,1)とも深く関わっています. 今,\(\small Z, W\) をN(0,1)にしたがう互いに独立な確率変数とします. それらの比 \(\small X=Z/W\) を考えると,\(\small X\) が したがう確率分布がコーシー分布 Cauchy(0,1) なのです. その理由は簡単です. \(\small W^2\) は自由度 1 のカイ自乗分布にしたがうので, \[\small X=\frac{Z}{\sqrt{W^2/1}}\] とみることができ,\(\small t\) 分布の定義より \(\small X\) は自由度1の \(\small t\) 分布にしたがいます. したがって,前述のことから,その確率密度関数は Cauchy(0,1) と一致します.
 以上のことを,N(0,1)からの乱数を利用して実際に確かめてみましょう. wxMaximaを利用します.

 (%i3)では,N(0,1)にしたがう乱数の比からなる1万個の要素を持つ リストを作成しています. (%i4)で,区間 \(\small [-10,10]\) を100等分したときの度数分布を 作成して,(%i6)でグラフ描画のときの書式を定めておいて, (%i7)によりグラフを描画しています.


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

平均と分散
 コーシー分布の平均と分散が存在しないことは,実際に積分してみれば分かります. 簡単のため, \(\small f(x,0,1)\) で考えます.

 (%i8)により,全範囲での積分の値は 1 になり, 確率密度関数の条件を満たしています.広義積分なので, 「ldefint」を利用しました.
 次に,\(\small xf(x,0,1), x^2f(x,0,1)\) の不定積分を求めると (%o11)(%o12)のようになり, いずれも \(\small x\to\pm\infty\) のときの極限値は存在しない ことが分かります.ここで,平均の計算式は \[\small \int_{-\infty}^{\infty}xf(x,0,1)\,dx =\int_{-\infty}^{\infty}\frac{x}{\pi(x^2+1)}\,dx\] で表されます.被積分関数が奇関数なので, 一見すると積分結果は 0 としがちですが, \[\small\begin{align*} (右辺) &=\lim_{\substack{b\to\infty\\ a\to-\infty}}\bigg[\frac1{2\pi}\log(x^2+1)\bigg]_a^b\\ &=\lim_{\substack{b\to\infty\\ a\to-\infty}}\bigg[\frac1{2\pi}\log\frac{b^2+1}{a^2+1}\bigg]_a^b \end{align*}\] となり,この積分の値(つまり平均)は存在しません. たとえば \(\small a=-tb~(t\gt 0)\) の場合は \[\small\begin{align*} (右辺) &=\lim_{b\to\infty}\frac1{2\pi}\log\frac{b^2+1}{t^2b^2+1}\\ &=\frac1{2\pi}\log\frac1{t^2}=-\frac1{\pi}\log{t} \end{align*}\] となるので,これは \(\small (-\infty, \infty)\) の全ての値を取り得ます. 平均が存在しないので,分散・歪度・尖度はいずれも存在しません.
 なお,(%i9)により確率密度関数の不定積分が存在するので, コーシー分布の累積分布関数は次のように逆正接関数を用いて 表されます. \[\begin{align*} \small {\rm cdf}\_{\rm cauchy}(x,0,1) &\small =\int_{-\infty}^{x}f(x,0,1)\,dx\\ &\small =\frac1{\pi}\big[{\rm arctan}(x)\big]_{-\infty}^{x}\\ &\small =\frac1{\pi}\left({\rm \arctan}(x)+\frac{\pi}{2}\right)\\ &\small =\frac12+\frac1{\pi}{\rm \arctan}(x) \end{align*}\]  したがって,この逆関数としての「quantile_cauchy(x,0,1)」は, \[\begin{gather*} \small \frac12+\frac1{\pi}{\rm \arctan}(y)=x\\ \small {\rm arctan}(y)=\pi\left(x-\frac{1}{2}\right)\\ \small \therefore\quad y=\tan\pi\left(x-\frac{1}{2}\right) \end{gather*}\] となり,正接関数になります. このことから,コーシー分布の四分位数の値は, たとえば次のようになります. \[\begin{align*} &\small {\rm quantile}\_{\rm cauchy}(1/4,0,1)=-1\\ &\small {\rm quantile}\_{\rm cauchy}(1/2,0,1)=0\\ &\small {\rm quantile}\_{\rm cauchy}(3/4,0,1)=1 \end{align*}\]

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

標本平均の分布
 このような分布にしたがう乱数の標本平均を取ると, それはどのようになっているのでしょうか.これまで見てきた分布では, 平均を取る個数を増やすと標本平均は母平均の周りに集まってきます. しかし,コーシー分布では母平均自体が存在しません. そこで,前と同じようにして,標本平均からなるリストを作って, その様子を調べてみましょう.

 (%i13)では,\(\small a=0, b=1\)の場合の コーシー分布にしたがう乱数のリストを作成しています. (%i14)では,そのリストの標本平均からなるリストを作成しています. コーシー分布のグラフから, 絶対値が 5 を越える確率は小さいはずなので, 標本平均の絶対値が 5 を越えるものの個数を数えてみましょう. (%i17)では,\(\small n\) 個の標本平均値からなるリスト hlist(n) の要素について, block関数を利用して 「\(\small i\) 番目の要素の絶対値が 5 を越えたら, kosu を 1 増やす」という操作を個々の要素ごとに調べています. 最初の [i, kosu] では,内部のローカル変数を指定します. そして,kosuに0を割り当てた後で,絶対値が5を越えたらkosuを1増やしす操作を 繰り返します.「for-while-do」では, 「,」で区切られる箇所までが「do」の内容です. 最後に,個数「kosu」を返してblock関数が終了します.
 このように,block関数を利用すると,繰り返しや条件分岐などを含む いろいろな処理を行わせることができます.

 (%i22)は,2個の平均からなる100要素のリストで,13要素の絶対値は 5を越えていることが分かります.(%i23)では標本平均を取る個数を10個に増やして5を 越える要素数を調べています. 標本平均を取る個数を増やしても,絶対値が 5 を越える要素数は同じになっています. ただし,乱数により計算しているので, たまたま個数が一致しただけです. 正規分布であれば,この個数は 「平均を取る個数を増やすと減る」はずです.

 ここで,幾つかの事項を,標準正規分布と比較してみましょう. 下図は,標準正規分布 N(0,1) のグラフとコーシー分布のグラフを重ね合わせています.

 次に,絶対値が 5 を越える確率はどの程度あるかを確かめてみましょう.

 (%i36)では,コーシー分布の分布関数「cdf_cauchy」を利用しています. ここでは,\(\small 2\times {\rm P} (X\leq -5)\) を求めています. 単純に1個抽出すると, \(\small {\rm P} (|X|\geq 5)=0.1256\cdots\) ということであり, けっこうな確率で大きな値が抽出されてしまいます. (%i37)は標準正規分布の場合です.この場合は, \(\small |X|\geq 5\) となることはほとんど無いといってよいでしょう.

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

 次に,標本平均を取る個数を増やすごとに,標本平均の値がどのように変化して いくかを調べてみましょう.

 (%i39)では,グラフを描くときの \(\small x\) 軸の値の リストを定めています.標本平均を取る個数を, 2 刻みで増やして 100 個の標本平均まで調べます. (%i40)は,コーシー分布にしたがう乱数の標本平均を, 2個の標本平均から100個の平均まで求めたリストです. (%i41)によりLyのグラフを描くと,下図のようになります. 単純に2個ずつ増やしながら標本平均を計算するとき, 標本平均の値の変化を表します.

 乱数の標本平均なのでグラフは描画するたびに異なりますが, 標本平均を取る個数を増やしいっても値が一定値に近づくような傾向は全く見られません. 下図は,10刻みで1000個の標本平均まで求めた場合です. 状況は,100個までの場合と全く変わりません. つまり,標本数を幾ら増やしても,その平均が幾らになるのか 「全く見当がつかない」ということです.おまけに, 左右対称の釣り鐘状の分布のはずなのに, 標本数を幾ら増やしても標本平均の絶対値が10を軽く超える場合が度々あるのです.
 同様のことを標準正規分布で調べるには, 「random_cauchy(0,1)」を「random_normal(0,1)」に 変更するだけです.自分で試してみてください.


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

 次に,作成した平均のヒストグラムをみてみましょう. 「TeXmacs+Maxima」ではヒストグラムの区分の際に エラーが生じるので,wxMaximaを利用します. コーシー分布にしたがう乱数「random_cauchy(0,1)」を利用して 下図の(%i16)により5個の乱数の平均からなる10000個の要素を持つリストを作成し, 区間 [-5, 5] を50等分してヒストグラムを描画したものです. 「draw」パッケージを利用して, (%i18)によりコーシー分布の確率密度曲線「pdf_cauchy(x,0,1)」も 同時に描画しました.

 この図を見ると,5個の標本平均の分布はもとのコーシー分布と 同じような曲線になっているようです. このことは一般に成り立ちます.つまり, コーシー分布にしたがう確率変数の標本平均は, 元のコーシー分布と同じ分布になるのです. 標本数を幾ら増やしても同じです. こちらも参照してください.
 以上のことから,コーシー分布は, 確率密度関数の形は正規分布と似ていますが, 正規分布とは全く異なる分布であることが分かります. コーシー分布のように, 平均も分散も存在しないような分布を含む確率分布に, 「安定分布」と呼ばれるものがあります. コーシー分布は安定分布です. 詳しくは「こちら」を参照してください.

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

参照したWebサイト
 コーシー分布の箇所では,下記のサイトを参照させていただきました.

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


copyright